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SUMMARY 


This  report  applies  the  method  of  canonical  eigenvalues  to  double-duct  propagation 
of  underwater  acoustics.  The  class  of  all  two-layer  bounded  profiles  with  free-  or  rigid¬ 
boundary  conditions  is  treated.  The  canonical  eigenvalue  equation  contains  two 
parameters.  (The  conventional  approach  has  six  parameters.)  Eigenvalue  solutions  are 
presented  for  23  sets  of  parameters  and  boundary  conditions.  These  solutions  can  be 
applied  to  any  equivalent  wave  problem  in  any  discipline  of  physics.  Physical  results 
are  presented  for  six  sound-speed  profiles.  These  include  plots  versus  frequency  of 
phase  velocity,  group  velocity,  and  eigenfunction  normalization  coefficients.  Eleven 
examples  of  eigenfunctions  (standing  waves)  are  presented.  At  the  critical  frequencies 
of  maximum  coupling  between  ducts,  the  phase  velocities  of  adjacent  modes  are 
pinched  toward  each  other,  the  group  velocity  and  eigenvalue  normalization  curves  for 
adjacent  modes  cross  each  other,  and  the  eigenfunction  amplitudes  of  adjacent  modes 
strongly  resemble  each  other.  The  report  develops  an  unexpected  example  in  which 
two  radically  different  profiles  have  identical  plots  of  phase  and  group  velocity  as  a 
function  of  frequency.  Special  configurations,  predicted  by  ray  theory  as  degenerate,  do 
not  lead  to  degenerate  eigenvalues.  The  theory  is  extended  to  treat  sound-speed  discon¬ 
tinuities  and  layers  with  different  constant  densities. 
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1.  INTRODUCTION 


This  report  applies  the  general  theory  of  canonical  eigenvalues  of  reference  1  to  the 
case  of  two-layer,  bounded,  underwater  acoustic  ducts.  The  general  canonical  eigen¬ 
value  equation  is  expressed  in  terms  of  two  mathematical  variables,  x  and  y,  and  a  set 
of  profile-related  parameters,  P,  and  Si+j.  Reference  2  treated  several  simple  sound- 
speed  profile  classes  for  which  there  are  only  one  mathematical  variable  (x  or  y)  and 
either  one  (Pi)  or  no  parameters.  Reference  3  treated  the  class  of  one-layer  bounded 
ducts  for  which  there  are  the  two  mathematical  variables,  x  and  y,  and  no  parameters. 

The  canonical  eigenvalue  equation  for  the  profile  classes  of  this  article  contains  the 
two  mathematical  variables  and  the  two  profile-related  parameters  Pj  and  S2.  These 
classes  of  profile  are  of  particular  interest  because  they  represent  the  simplest  configu¬ 
ration,  which  exhibits  the  coupling  characteristics  of  double-duct  propagation.  Our  first 
extensive  treatment  of  double  ducts  is  presented  in  reference  4,  which  involved  an 
eigenvalue  matrix  of  rank  nine.  The  canonical  eigenvalue  approach  was  developed  to 
simplify  the  analysis  of  double-duct  propagation.  Thus,  the  present  report  represents 
our  first  application  of  the  method  to  the  problem  for  which  it  was  developed. 

Appendix  B  extends  the  general  theory  of  reference  1  to  include  the  treatment  of 
sound-speed  discontinuities  at  layer  interfaces  and  the  treatment  of  layers  of  different 
constant  density.  Specific  computing  forms  are  given  for  the  two-layer  profile  bounded 
above  by  a  free  surface  and  below  by  a  rigid  surface.  However,  this  theory  has  not 
been  programmed  as  it  is  beyond  the  scope  of  the  present  effort.  The  number  of 
canonical  parameters  increases  to  four,  which  is  twice  the  number  for  continuous  pro¬ 
files  and  no  treatment  of  density. 
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2.  THEORY 


This  major  section  and  the  appendices  contain  most  of  the  theoretical  treatment  of 
this  report.  The  equations,  which  follow  directly  from  the  application  of  reference  1  to 
the  profile  class  of  figure  1,  represent  a  minor  part  of  the  treatment.  The  equations 
necessary  to  generate  the  numerical  examples  of  this  report  are  presented  as  well  as 
the  analysis  thereof.  Some  of  this  theory  may  appear  arcane  on  first  exposure.  How¬ 
ever,  it  will  become  clearer  when  examined  later  in  the  context  of  the  relevant  numeri¬ 
cal  examples  as  presented  in  section  3. 

A.  PROFILE  REPRESENTATION 

The  sound-speed  profile  consists  of  two  layers.  The  sound  speed  in  each  of  these 
layers  is  expressed  as 


[Ci/C(Z)]2  =  l-27i(Z-Zi)/Ci  (1) 

A  TWO-LAYER,  BOUNDED,  CANONICAL  PROFILE 

CASE  A  CASE  B 


Figure  1.  Schematic  of  the  profile  class. 

where  Q,  Zj,  and  7]  are  the  sound  speed,  depth,  and  sound-speed  gradient,  respec¬ 
tively  at  the  top  of  layer  i. 

Figure  1  presents  a  schematic  of  the  profile  class.  The  profile  class  is  characterized 
by  two  parameters,  Pj  and  S2,  related  to  the  parameters  of  eq.  1  by 

Pi  =  (~710/72)1/3  (2) 

and 

S2  =  [(C1/C3)2  -  l]/[(Cj/C2)2  -  1]  .  (3) 

Here,  Q  is  the  surface  sound  speed,  C2  is  the  sound  speed  at  the  interface  between 
the  two  layers,  and  C3  is  the  bottom  sound  speed.  In  eq.  2,  72  is  the  gradient  at  the 
top  of  layer  2;  whereas  7io  is  the  gradient  at  the  bottom  of  layer  1  as  given  by 

^10  =  C23  Tj/C,3  (4) 

where  7i  is  the  gradient  at  the  surface. 
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As  demonstrated  in  reference  1,  the  parameterization  of  eqs.  2  and  3  is  uniquely 
determined  by  the  condition  that  the  frequency  and  the  phase  velocity  do  not  appear  in 
the  canonical  eigenvalue  equation. 

The  sound-speed  profile  is  completely  determined  by  give  profile  parameters.  These 
are  Ci,  C2,  C3,  7i,  and  72.  However,  the  customary  eigenvalue  equation  usually 
replaces  72  by  the  dimensionless  parameter  Pj.  Here,  12  can  be  related  to  Pi,  Ci,  C2, 
and  7i  through  eqs.  2  and  4. 

In  the  canonical  eigenvalue  approach,  the  eigenvalue  solutions  are  obtained  for  a 
fixed  pair  of  values  of  the  parameters  p  and  S.  (From  this  point  on  we  will  drop  the 
subscripts  on  Pi  and  S2  to  simplify  the  notation).  Each  point  in  the  p,  S  plane  repre¬ 
sents  a  different  profile  class.  To  translate  the  canonical  eigenvalues  from 
mathematical  space  to  physical  space,  we  must  specify  a  particular  sound-speed  profile 
of  the  class.  Here,  we  are  at  liberty  to  select  Cj,  C2,  and  7).  The  parameter  C3  may 
then  be  determined  from 

C3  =  C1{[(CI/C2)2-l]S  +  ir1/2  .  (5) 

Thus,  the  profile  is  completely  specified  by  the  three  selected  parameters,  eq.  5,  and  p. 

Cases  A  and  case  B  are  distinguished  in  figure  1,  because  the  canonical  eigenvalue 
equation  yields  two  sets  of  modes.  The  set  of  modes  corresponding  to  case  A  are  of 
most  interest  because  they  exhibit  the  coupling  between  double  ducts.  For  case  A,  7i 
is  positive  and  C2  is  greater  than  Q.  The  set  of  modes  corresponding  to  case  B  does 
not  exhibit  any  coupling  phenomenon.  For  case  B,  7i  is  negative  and  C2  is  smaller 
than  Ci. 

Next  we  examine  p  and  S  in  the  context  of  figure  1.  Our  interest  is  this  article  is 
limited  to  the  base  of  p  positive.  This  means  that  the  gradient  of  the  two  layers  has 
opposite  signs.  If  p  is  negative,  then  there  is  no  double-duct  configuration,  because  C3 
of  case  A  is  larger  than  C2. 

Next,  consider  the  allowable  -  nge  of  values  on  S  for  which  there  are  double  ducts 
in  case  A.  When  Q  =  C3,  S  =  0.  When  C2  =  C3,  5  =  1.  When  C3  =  0,  S2  =  -  00.  Our 
primary  interest  is  in  0  <  S  <  1  for  which  Q<  C3  <  C2  and  in  modestly  negative 
values  of  S  for  which  C3  is  somewhat  less  than  Ci. 

We  examined  alternative  representations  for  S  in  terms  of  Zj  and  Z^  where  Z]  and 
Z2  is  the  thickness  of  layers  1  and  2  respectively.  We  obtained 


S  =  (C2/C3)2  (1-p-3Z2/Z,)  . 


(6) 
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This  does  not  appear  a  particularly  useful,  with  the  exception  of  the  case  where  S  =  0. 
Here,  eq.  6  reduces  to 

Z2=p3Z!.  (7) 

Our  main  interest  is  in  the  free  surface,  rigid-bottom  boundary  condition  in  fig¬ 
ure  1.  However,  we  shall  also  consider  free-free  and  rigid-rigid  boundary  conditions  as 
they  illustrate  certain  features  of  interest  in  the  solution.  (We  need  not  consider  the 
rigid-free  boundary  condition  because  this  can  be  obtained  by  inverting  the  free-rigid 
condition.)  We  note  that  the  analysis  of  the  sound-speed  profile  presented  here  applies 
equally  well  to  all  boundary  conditions.  The  boundary  conditions  determine  the  particu¬ 
lar  form  of  the  eigenvalue  equation,  but  have  nothing  to  do  with  the  parameters  p  and 
5  and  their  relation  to  the  profile  parameters. 


B.  CANONICAL  EIGENVALUE  FORMULATION 

First,  consider  the  eigenvalue  formulation  for  the  free-rigid  boundary  condition.  The 
eigenvalue  equation  may  be  expressed  as 


Ai(-  £n) 
Ai(-  £io) 

P  Ai'(-  £io) 
0 


Bi(-Sn) 
Bi(-  S,0) 
PBi'(-U) 
0 


0 

-  Ai(-p2£10) 

Ai'(-P2£io) 

Ai'(-  £20) 


0 

-  Bi(-p2£]0) 

Bi'(~P2U) 

Bi'(-ko) 


=  0  .  (8) 


Rows  1  to  4  of  the  matrix  follow  respectively  from  eqs.  23,  18,  22,  and  24  of  refer¬ 
ence  2. 

We  now  consider  the  transformations  into  the  canonical  variables  of  reference  1 
given  by 


£11  =  x 


(9) 


and 


£io  =  y-  0°) 

One  of  the  canonical  parameters  is  p.  The  other,  S,  enters  the  eigenvalue  equation 
through 

£20  =P2(1  -  S)x  +p2Sy  (11) 

which  follows  from  eqs.  32  and  37  of  reference  1. 
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Next,  we  introduce  various  sets  of  equations  that  facilitate  the  expansion  of  the 
matrix  of  eq.  8  and  its  solution  by  Newton’s  method.  Certain  combinations  of  Airy 
functions  associated  withe  the  interface  between  layers  1  and  2  are  common  to  the 
solution  for  any  combination  of  boundary  conditions.  We  define  these  combinations  as 
Gi  to  G4  where 


G!  =  Ai(-  y)Ai'(-p2y)  +  pAi'(-  y)Ai(-p2y) 

(12) 

G2  =  Bi(-  y)  Ai’  (-  p2y)  +  pBi'  (-  y)  Ai(-  p2y) 

(13) 

G3  =  Ai(-  y)Bi'(-p2y)  +  pAi'(-  y)Bi(-p2y) 

I 

(14) 

G4  =  Bi(-  y)Bi'(-p2y)  +  pBi'(- y)Bi(-p2y)  . 

(15) 

We  shall  also  need  the  derivative  of  these  functions  with  respect  to  y,  holding  p 
fixed.  These  derivations  are 

Gi*  =  -  (1  +  p>)[Ai'(-y)Ai'(-p2y)  -py  Ai(-  y)Ai(-  f?y)]  (16) 

G2'  =  -  (1  +p3) [Bi'(-  y)Ai'(-p2y)  -py  Bi(- y)Ai(-p2y)]  (17) 

G3'  =  -  (1  +P3)[Ai'(-y)Bi'(-p2y)  -py  Ai(- y)Bi(-p2y)]  (18) 

and 

G4'  =  -  (1  +p3)[Bi'(-y)Bi'(-p2y)  -py  Bi(-  y)Bi(-p2y)]  .  (19) 

We  note  in  passing  the  Gi  and  Gi'  appear  in  eqs.  79  and  84  of  reference  2  as  the  so¬ 
lution  to  the  unbounded  refractive  duct.  Moreover,  G2  appears  in  eq.  105  of  reference 
2  as  part  of  the  the  canonical  eigenvalue  equation  for  a  refractive  duct  bounded  above 
by  a  free  surface. 

Next,  we  define  certain  combinations  of  G]  to  G4  as  given  by 


Gn  =  Bi'(-  £20)  G2  -  Ai'(-  C20)  G4 

(20) 

G12  =  Bi'(-£20)  Gj  -  Ai'(-C2o)G3 

(21) 

G,3  =  Bi(-  ^2o)  G2  -  Ai(-  ^2o)  G4 

(22) 

and 

G14  =  Bi(-C2o)  G2  -  Ai(-U)  G4. 

(23) 
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We  note  that  £20  is  given  in  terms  of  x  and  y  by  eq.  11.  However,  we  retain  the  use  of 
£20  here  to  simplify  the  manner  of  writing  the  expressions. 

Next,  we  define  certain  combination  of  Gn  to  Gi4  as  given  by 


G15  =  Ai(— x)  Gn  -  Bi(-x)Gj2 

(24) 

G16  =  Ai'(-x)G„  -  Bi'(-x)  G32 

(25) 

G17  =  Ai(-x)  G!3  -  Bi(-  x)  Gi4 

1 

(26) 

G18  =  Ai'(-  x)  G13  -  Bi'(-x)  G14. 

(27) 

We  also  require  the  quantities  G15',  G^',  and  Gn',  which  are  defined  respectively  as 
Gis,  Gi6,  and  Gn  with  Gi,  G2,  G3,  and  G4 replaced  respectively  by  Gj',  G2',  G3',  and 
G4'. 

We  are  now  ready  to  present  the  eigenvalue  equation  and  other  expressions  used  in 
the  solution.  The  canonical  eigenvalue  equation  for  the  free-rigid  boundary  is  given  by 

Gis  =  0  .  (28) 

Equation  28  can  be  obtained  by  expanding  the  matrix  of  eq.  8  by  the  elements  of  the 
first  row.  Equation  28  is  solved  by  Newton’s  method  similar  to  that  used  in  reference  3 
for  one-layer  bounded  ducts.  This  method  requires  one  of  the  two  derivatives 

<5  Gls/<5  x  =  -  G,e  +  WO  -  S)  G17  (29) 


or 


<5  Gi5/<3  y  =  -  G'is  +  £20p2  S  G]7  .  (30) 

Item  Gi6  arises  from  the  Airy  function  with  a  -x  argument.  Item  G']5  arises  from  the 
Airy  functions  with  a  -y  or  -p*y  argument.  Item  GJ7  arises  from  the  Airy  function  with 
a  -£>(,  argument.  One  solves  eq.  28  by  either  selecting  a  fixed  value  of  y  and  iterating 
to  determine  x  by  eq.  29  or  by  selecting  a  fixed  value  of  x  and  iterating  to  determine 
y  by  eq.  30. 

Next,  consider  the  free-free  boundary  condition.  For  this  case,  the  Airy  function 
derivatives  of  row  4  of  eq.  8  are  replaced  by  the  Airy  functions.  The  eigenvalue  equa¬ 
tion  is  given  by 

Gn  =  0  .  (31) 
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with  derivatives  given  by 


<5G,7/<5x  =  -G18-p2(1-S)  G,s 


(32) 


and 


5G„/dy  =  -G„’  -p2SG,5. 


(33) 


In  the  case  of  the  rigid-rigid  boundary  condition,  the  Airy  functions  of  row  1  of 
eq.  8  are  replaced  by  the  Airy  function  derivatives.  The  eigenvalue  equation  is  given 
by 


Gi6  =  0  . 


(34) 


with  derivatives  given  by 


<5G16/<5x  =  -G15  +  W(1-S)  G18 


(35) 


and 


<5  G16/<5  y  =  -  G')6  +  ?2oP2  S  Gig  . 


(36) 


We  note  that  it  is  relatively  easy  to  include  the  other  boundary  conditions  in  the 
analysis  for  Gi?  and  Gi6,  the  eigenvalue  expressions  for  the  free-free  and  rigid-rigid 
conditions  are  also  required  in  eq.  29  for  use  in  solving  eq.  28  for  the  free-rigid  condi¬ 
tion.  Our  computer  implementation  consists  of  implementing  eqs.  12  to  27  in  a  general 
program  and  using  them  as  needed  in  eqs.  28  to  30  or  their  other  boundary  condition 
counterparts. 

C.  ASYMPTOTIC  SOLUTIONS 

This  section  is  concerned  with  the  eigenvalue  solutions  for  infinite  negative  y  or 
infinite  negative  x. 

First,  we  consider  the  behavior  for  large  negative  values  of  y.  For  all  three  bound¬ 
ary  conditions,  the  dominant  term  is  G4of  eq.  15  because  it  contains  products  of  Airy 
functions,  both  factors  of  which  become  infinite.  If  we  now  divide  the  eigenvalue  equa¬ 
tion  through  by  G4,  all  of  the  terms  will  go  to  zero  except  the  one  that  contained  G4 
before  the  division.  It  follows  then,  that  the  eigenvalue  can  be  reduced  to 


Ai(-  x)  Ai'(-  £20)  =  0 

(37) 

Ai(-  x)  Ai(-  ?20)  =  0 

(38) 

8 


and 


Ai'(-x)  Ai'(-  £2o)  =  0 


(39) 


for  the  free-rigid,  free-free,  and  rigid-rigid  conditions  respectively. 

Let  the  negatives  of  the  zeroes  of  Ai  and  Ai'  be  given  by  an  and  a'n.  These  zeroes 
are  tabulated  in  reference  5.  In  the  case  of  a  free  surface,  there  are  vertical  asymp¬ 
totes  in  an  x,y  plot  at  y  =  -  00  for 

x  =  an  .  (40) 

In  the  case  of  a  rigid  surface,  there  are  vertical  asymptotes  at  y  =  -  °°  for 

x  =  a'„  .  (41) 

In  the  case  of  a  rigid  bottom,  the  variable  y  is  asymptotic  to  the  lines 

y  =  (1  -S_1)x  +  a'n/p2  S,  S  *  0  .  (42) 

In  the  case  of  S  =  0,  there  are  vertical  asymptotes  for 


x  =  a'n/ p2 


(43) 


In  the  case  of  a  free  bottom,  the  counterparts  of  eqs.  42  and  43  are 
y  =  (1  -  S_1)x  +  an/ p2  S,  S  *  0 


(44) 


and 


x  =  an/p2  .  (45) 

As  will  be  demonstrated  in  section  3,  the  curves  with  negative  values  of  y  corre¬ 
spond  to  the  case  A  modes  of  figure  1.  The  asymptotic  solutions  of  eqs.  40  to  45  are 
of  particular  importance  because  they  provide  a  first  approximation  to  the  critical 
frequencies  of  double-duct  propagation.  The  asymptotic  solutions  are  also  important  as 
a  starting  point  for  the  iterative  solutions  of  the  eigenvalue  equation. 

Next,  consider  the  behavior  for  large  negative  values  of  x.  In  the  case  of  eq.  28, 
the  dominant  term  is  -Bi(-x)  Bi'(-^0)  G,.  Here,  the  first  two  factors  both  become 
infinite  when  S<1.  If  we  divide  the  eigenvalue  equation  by  these  factors,  all  terms  go 
to  zero  with  the  exception  of  G,.  Thus,  here  there  are  horizontal  asymptotes  at 
x  =  -  00  for  the  values  of  y  which  satisfy 
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G1  =  0  . 


(46) 


Table  3  and  figure  6  of  reference  2  gives  the  solutions  of  y  as  a  function  of  p.  If  we 
examine  the  other  boundary  conditions,  vve  find  that  the  dominant  term  is  again  G, 
multiplied  by  a  product  of  Bi  or  Bi'.  The  result  is  that  eq.  46  applies  to  all  of  the 
boundary  conditions. 

We  note  that  x  equal  to  minus  infinity  corresponds  to  infinite  frequency  for  the 
case  B  modes  of  figure  1.  We  also  note  that  eq.  46  is  identical  to  the  solution  for  an 
unbounded  refractive  duct  as  given  by  eq.  79  of  reference  2.  At  sufficiently  high  fre¬ 
quencies,  the  solution  for  case  B  of  figure  1  is  independent  of  the  boundary  conditions 
and  approaches  the  solution  for  the  unbounded  refractive  duct  as  a  limit. 

D.  LATTICE  POINTS 

In  the  case  of  the  ducts  of  reference  3,  certain  finite  roots  of  the  eigenvalue  equa¬ 
tion  could  be  related  to  the  zeros  of  the  various  Airy  functions.  These  zeros  formed 
sets  of  lattice  points  through  which  the  eigenvalue  curves  pass.  Next,  we  examine  the 
possibility  for  our  two-layer  ducts.  First,  consider  the  case  of  S  =  0  and  p  =  1.  For 
these  values,  eq.  31,  for  the  free-free  boundary,  is  satisfied  if  both  Ai(-x)  and 
G,(y)  =  0.  Thus,  here  the  eigenvalue  curves  go  through  the  lattice  points  given  by 

x  =  an;  y  =  am  or  y  =  a'm  .  (47) 

This  notation  means  that  the  zero  Ai(-y)  or  Ai'(-y)  can  be  of  a  different  order  than 
that  of  Ai(-x). 

We  note  that  Ai(-x)  and  Ai(-y)  or  Ai'(-y)  also  go  to  zero  at  x  or  y  =  -  ».  Thus, 
these  points  at  infinity  are  also  lattice  points  of  eq.  47.  They  take  the  form  of  horizon¬ 
tal  asymptotes  for  x  -*  -  °°  and  vertical  asymptotes  for  y  -*  -  oo . 

Our  previous  analysis  leading  to  eq.  46,  resulted  in  the  same  conclusion.  The  dif¬ 
ference  is  that  eq.  46  applies  to  all  cases,  whereas  eq.  47  applies  only  to  the  case  of 
S  =  0  and  p  =  1  for  the  free-free  boundary  condition. 

We  find  that  eq.  31  is  also  satisfied  if  both  Bi(-x)  =  0  and  G4  =  0.  Here,  the 
eigenvalue  curves  go  through  the  lattice  points  given  by 

x  =  bn;  y  =  bm  or  y  =  b'm  .  (48) 

where  bn  and  b'n  represent  the  negative  of  zeros  of  Bi  and  Bi'.  We  note  that  this  set 
of  lattice  points  does  not  contain  the  points  at  x  or  y  =  -  00  because  Bi  or  Bi'  do  not 
go  to  zero. 


10 


Next,  consider  the  case  of  the  rigid-rigid  boundary  conditions  as  applies  to  eq.  34 
for  S  =  0  and  p  =  1.  The  counterparts  of  eqs.  47  and  48  are 

x  =  a'n;  y  =  am  or  y  =  a'm  .  (49) 

and 

x  =  b'„;  y  =  bm  or  y  =  b'm  .  (50) 

Consider  now  what  happens  if  we  set  S  =  0  and  p  =  1  in  eq.  28  for  the  free-rigid 
boundary  condition.  The  combination  of  functions  corresponding  to  eq.  47  is 

x  =  an;  y  =  yn  (51) 

where  yn  is  a  solution  of 

Bi'(-a„)G,(y)  -  Ai'(-an)  G3(y)  .  (52) 

The  problem  here  is  that  the  zeros  depend  on  an .  The  solution  of  eq.  52  only  gives  us 
one  point  in  the  x,y  plane  rather  than  the  lattice  of  points  given  by  eqs.  47  to  50.  The 
same  thing  happens  if  we  take  S  *  0  or  p  *  1  for  any  of  the  boundary  conditions.  We 
can  simplify  the  equation  somewhat  for  x  =  an,  bn,  a'n,  or  b'n.  However,  this  leads  to 
an  equation  in  y  which  contains  functions  of  an,  bn,  a'n,  or  b'n  and  there  is  no  signifi¬ 
cant  simplification. 

For  the  free-free  and  rigid-rigid  ducts  of  reference  3,  the  line  y  =  x  was  a  solution 
to  the  eigenvalue  equation.  If  we  substitute  y  =  x,  £20  reduces  to  p*x  and  G17  of  eq.  26 
and  Gie  of  eq.  25  both  identically  go  to  zero.  As  we  shall  see  later,  this  solution  plays 
a  role  in  the  eigenvalue  plots. 

E.  EXPRESSIONS  FOR  FREQUENCY  AND  PHASE  AND  GROUP 
VELOCITY 

It  is  convenient  now  to  introduce  expressions  that  convert  the  canonical  variables 
into  physical  quantities  of  interest.  As  discussed  in  the  derivation  of  eq.  14  of  refer¬ 
ence  1,  the  frequency  may  be  written  as 

f  =  K  fca  (53) 


where 

fca=|y-x|3/2 


(54) 
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and 


K  =  |C]2  C2'2  -  1|"3/2  lyiK'1  •  (55) 

Here,  fca  of  eq.  54  is  known  as  the  canonical  frequency.  It  is  related  to  the  physical 
frequency  by  the  scale  constant  K  of  eq.  55.  As  expressed  in  reference  1  for  the 
case  A  profile  y-x  and  C12C2"2  -  1  are  both  positive  while  both  are  negative  for  the 
case  B  pufile.  This  is  how  the  two  sets  of  eigenvalues  are  related  to  their  respective 
profiles.  In  any  case,  the  use  of  absolute  value  signs  in  eqs.  54  and  55,  allows  us  to 
treat  both  cases  A  and  B  modes. 

Equation  53  gives  the  frequency  as  a  function  of  x,  since  the  eigenvalue  solutions 
gives  y  as  a  function  of  x.  The  phase  velocity  may  now  be  determined  from 

Cp  =  C,[l  -  f-2/37r2/3|yi|2/3x]-i/2  (56) 

This  equation  follows  from  the  discussion  of  eq.  17  of  reference  1.  Equation  56  gives 
the  phase  velocity  as  a  function  of  frequency. 

The  expression  for  group  velocity  is  given  by 

Cg  =  Cp[l  +  (Cp2  -  C]2  )/3Cj2  +  Cp2  (C2~2  -C,"2  )L/3]-J  (57) 


where 

L  =  (1  -  dy/dx)'1  .  (58) 

The  deviation  of  eq.  57  is  given  in  the  discussion  of  eq.  69  in  reference  1.  Here,  dy/dx 
is  the  slope  of  the  eigenvalue  curves  and  is  evaluated  as 

dy/dx  =  -  (<5G/<5x)/(<5G/<5y)  (59) 

with  the  use  of  partial  derivatives  of  section  2B. 

Equations  53  to  59  can  be  applied  to  all  canonical  eigenvalue  problems.  Our  appli¬ 
cation  to  two-layer  bounded  ducts  only  enters  the  equation  through  x  and  y  and  the 
particular  form  used  in  eq.  59. 

F.  ESTIMATE  OF  CRITICAL  FREQUENCIES 

An  important  aspect  of  the  coupling  between  double  ducts  are  the  critical  frequen¬ 
cies.  These  are  the  frequencies  where  the  difference  in  phase  velocities  between 
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adjacent  modes  form  a  relative  minimum.  At  these  frequencies,  the  eigenfunctions  for 
adjacent  modes  display  a  similarity  associated  with  strong  coupling  between  the  ducts. 

Exact  evaluation  of  the  critical  frequencies  requires  a  detailed  evaluation  of  the 
mode  eigenvalues.  However,  certain  relatively  simple  approximations  to  these  critical 
frequencies  are  available  from  conventional  eigenvalue  approaches.  This  section 
reviews  the  results  of  the  conventional  approach,  develops  approximations  for  the 
canonical  eigenvalue  approach,  and  compares  the  results  for  the  two  approaches. 

In  the  case  of  conventional  mode  theory,  we  consider  plots  of  phase  velocity  versus 
frequency  for  each  of  the  idealized  single  ducts  of  the  double-duct  profile.  The  fre¬ 
quencies  where  these  plots  intersect  are  an  estimate  of  the  critical  frequencies  of  the 
double  duct.  For  over  20  years  we  have  used  the  phase-integral  method  of  reference  2 
to  obtain  these  estimates.  However,  the  work  of  reference  2  has  demonstrated  that  the 
exact  normal  mode  sc'utions  for  the  idealized  single  ducts  lead  to  better  estimates  than 
does  the  phase-integral  method. 

Consider  this  approach  for  the  free-rigid  boundary  condition  in  figure  1.  From 
eq.  54  of  reference  2,  the  phase  velocities  in  the  upper  duct  may  be  written  as 


(60) 


and  those  of  the  lower  duct  as 

Cp  =  C3[l  -  f-2/31T-2/3|y20|2/3a'in]-V2  (61) 

where  Y20  is  the  gradient  at  the  bottom  interface.  If  we  equate  eq.  60  to  61  and  solve 
for  f,  we  obtain 

f  =  Tr",|y)|(p'2a'm-an)3/2  (C^Cf2  - 1)3/2.  (62) 

To  obtain  eq.  62  we  make  use  of 

■  (63) 

Equation  63  follows  from  eq.  5  of  reference  1.  For  f  to  be  real 

P~2  a'm  >  an  for  Q  >  C3  (64) 


and 


P~2  am'  <  an  for  C,  <  C3  . 


(65) 
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Here,  m  and  n  are  independent  and  may  be  any  mode  numbers  provided  eqs.  64  or  65 
are  satisfied. 

Note  that  eq.  62  breaks  down  when  C,  =  C3.  This  means  that  the  curves  of  eqs.  60 
and  61  do  not  intersect  but  run  parallel  to  each  other  in  a  manner  similar  to  the 
curves  for  different  mode  numbers  in  a  single  duct.  When  C,  =  C3,  two  cases  can 
occur.  When 

|y,l2/3  a„  *  |y20|2'3  a'm  (67) 

i.e.,  when 


P2  =  a'm/an  ,  (68) 

the  curves  for  modes  m  and  n  are  superimposed. 

Equations  60  to  68  apply  to  the  free-rigid  boundary  condition.  However,  they  also 
apply  to  the  free-free  or  rigid-rigid  boundary  conditions  if  we  replace  a'm  by  am  or  an 
by  a'n,  respectively. 

Two  subcases  of  eq.  68  occur.  One  subcase  is  p  =  1  for  the  free-free  or  rigid-rigid 
boundary  condition.  Here,  the  counterparts  of  eq.  68  are  satisfied  for  m  =  n.  This  is  a 
case  of  symmetry  in  which  the  curves  of  eq.  60  are  identical  to  those  of  eq.  61.  Each 
curve  of  eq.  60  lies  on  the  corresponding  curve  of  eq.  61.  This  subcase  cannot  occur 
for  the  free-rigid  boundary  condition  because  the  mixed  boundary  condition  breaks  the 
profile  symmetry. 

The  second  subcase  occurs  when  there  is  no  profile  symmetry  and  the  gradients  are 
constructed  to  satisfy  eq.  68  or  its  free-free  or  rigid-rigid  counterparts  for  a  particular 
choice  of  m  and  n.  Here,  mode  n  of  the  upper  duct  is  superimposed  on  mode  m  of 
the  lower  duct.  In  the  case  of  the  free-rigid  boundary  condition,  m  may  be  equal  to  n. 

Equations  62  and  68  apply  to  the  double-duct  configuration  of  figure  1.  However, 
these  equations  have  similar  counterparts  for  any  double-duct  configuration.  There  is 
an  equivalent  to  eq.  62  so  long  as  the  minimum  sound  speed  for  each  of  the  two  ducts 
is  not  the  same.  Here,  critical  frequencies  are  formed.  An  example  of  this  case  is 
presented  in  reference  2  with  more  detail  presented  in  reference  4.  Here,  the  double¬ 
duct  configuration  is  formed  by  two  refractive  ducts,  each  with  its  own  axial  sound 
speed.  There  are  also  counterpart  equations  to  eqs.  66  and  68  when  the  minimum 
sound  speed  of  the  two  ducts  is  the  same. 

An  example  of  the  situation,  illustrated  by  eq.  68,  has  been  investigated  by  conven¬ 
tional  mode  approaches  in  an  unpublished  manuscript.  Some  limited  results  on  the 
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sound-speed  profile  and  phase  and  group  velocity  are  published  in  reference  6.  This 
investigation  was  carried  out  because  we  anticipated  that  superimposed  phase  velocity 
curves  would  represent  a  case  of  “supercoupling,”  since  strong  coupling  effects  were 
present  in  the  usual  case  where  the  phase  velocity  curves  intersected  at  a  single  point. 

The  profile  configuration  chosen  for  this  investigation  was  a  positive-gradient  sur¬ 
face  duct  overlaying  a  symmetric  refractive  duct.  Here,  eq.  60  applies  as  is.  Equa¬ 
tion  61  is  modified  to  accommodate  a  symmetric  refractive  duct  with  axial  gradients  of 
±  ya.  In  eq.  64,  C3  the  axial  sound  speed  is  equal  to  C,  and  ^20  is  replaced  by  ya. 
From  eq.  82  of  reference  2,  we  find  that  that  a'm  applies  for  the  odd  modes,  but  am 
applies  for  the  even  modes.  We  wished  to  superimpose  mode  1  of  the  upper  duct  on 
mode  1  of  the  lower  duct.  The  equivalent  to  eq.  68  then  becomes 

Y=Y  {a'\/a\?n  =  Y  x  0.2876287023  .  (69) 

1  «l  <1 

The  numerical  results  obtained  for  this  profile  will  be  discussed  in  section  3  when  an 
example  of  the  application  of  eq.  68  is  presented. 

Now,  consider  the  canonical  eigenvalue  approach.  As  will  .be  demonstrated  later, 
the  canonical  eigenvalue  plots  indicate  that  the  critical  frequencies  are  approximated  by 
the  intersection  of  the  vertical  asymptotes  of  eqs.  40  or  41  with  the  slant  asymptotes  of 
eqs.  42  or  44.  In  the  case  of  the  free-rigid  boundary  condition  at  the  point  of  intersec¬ 
tion,  eqs.  40  and  42  lead  to 

y  -  x  =  S_1(p'2  a'm  -  an)  .  (70) 

Our  interest  is  in  the  case  A  modes  for  which  y  <  x  and  C^  <  C2.  The  condition  that 
y  <  x  (corresponding  to  real  frequencies)  leads  to  the  same  conditions  as  eqs.  64 
and  65.  If  we  substitute  eq.  70  into  eq.  53,  we  obtain  eq.  62.  Thus,  the  intersection  of 
asymptotes  in  the  canonical  approach  leads  to  the  same  approximations  for  the  critical 
frequencies  obtained  from  customary  mode  theory  by  the  intersection  of  the  phase 
velocity  versus  frequency  curves  for  the  two  single  ducts. 

Next,  consider  the  case  of  S  =  0.  Here,  the  two  asymptotes  are  vertical  as  given  by 
eqs.  40  and  43.  These  asymptotes  are  superimposed  when 

p2=a'm/an,  (71) 

which  is  the  same  condition  as  eq.  68.  For  the  free-free,  or  rigid-rigid  cases,  eqs.  70 
and  71  apply  with  the  replacement  of  a'm  by  am,  respectively.  For  these  cases,  eq.  71 
is  satisfied  for  m  =  n  and  p  =  1,  which  is  the  same  circumstances  as  discussed  in  con¬ 
nection  with  eq.  68  in  the  conventional  approach. 
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Thus,  the  canonical  approximation  gives  the  same  result  as  the  conventional 
approximation  for  conditions  under  which  the  mode  eigenvalues  are  separated  or  not 
separated  for  the  case  of  S  =  0.  However,  there  is  a  significant  difference  between  the 
two  approximations  when  we  consider  frequency.  As  will  be  demonstrated  later,  two 
canonical  eigenvalue  curves  have  the  common  asymptote  of  eq.  71  or  its  equivalent  for 
the  other  boundary  conditions.  These  two  curves  only  intersect  at  y  =  -  oo  and  finite  x. 
From  eqs.  53  and  54,  this  corresponds  to  infinite  frequency.  Thus,  the  canonical 
approximation  yields  the  same  phase  velocity  for  these  paired  modes  only  at  infinite 
frequency;  whereas  the  conventional  approximation  states  that  the  paired  modes  have 
the  same  phase  velocity  for  all  frequencies. 

G.  EIGENFUNCTIONS 

This  section  treats  the  eigenfunctions.  Eigenfunctions  are  defined  here  as  solutions 
that  satisfy  the  Helmholtz  equation  and  the  boundary  and  interface  conditions. 

The  eigenfunction  in  layer  i  is  given  by 


Fi(Z)  =  Dj  Aj(-  £0+  Ej  Bi(-  £j)  .  (72) 

Equation  72  is  the  general  solution  to  the  Helmholtz  equation  which  for  the  profile  of 
eq.  1  reduces  to  Airy’s  differential  equation.  Here 


Ci  =  x-a,  Z 

(73) 

£2  =  P2y  -  p  1  aj(Z  -  Z2) 

(74) 

a,3  =-2yi<u2/Ci3  • 

(75) 

The  remaining  task  is  to  specify  Dj  and  E;  so  that  the  boundary  and  interface  con¬ 
ditions  are  satisfied.  First,  consider  the  case  of  a  free  surface.  Here, 


Di  =  Bi(-  x)  and  Ei  =  -  Ai(-  x)  (76) 

satisfy  the  surface  boundary  condition.  We  now  define  Fi(-y)  as  eq.  72  evaluated  at 
the  bottom  of  layer  1. 

F i (—  y)  =  Bi(-  x)  Ai(-  y)  -  Ai(-  x)  Bi(-  y)  .  (77) 

Next,  we  define 

F,'=dF/D£.  (78) 
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Hence, 


F,'(- y)  =  -  Bi(- x)  Ai'(- y)  +  Ai(- x)  Bi'(-  y)  .  (79) 

Equations  79  and  80  of  reference  1  give  the  solutions  for  D2  and  E2  as 
D2  =  ir  [Fj(-  y)  Bj'(-  p2y)  -pFi'(-  y)  Bi(-  p2y)]  .  (80) 

and 

Ei  =  -it  [Fj(-y)  A/(-  p2y)  -  pFi'(-  y)  Ai(-p2y)]  .  (81) 

Equations  80  and  81  satisfy  the  continuity  conditions  at  the  interface  between  layers  1 
and  2.  The  boundary  condition  at  the  bottom  is  satisfied  by  the  eigenvalue  equation, 
since  eqs.  76,  80,  and  81  only  ^atisfy  the  surface  and  interface  conditions. 

Now,  consider  the  case  of  a  rigid  surface.  Here,  we  replace  Ai(-x)  and  Bi(-x)  in 
eqs.  76,  77,  and  79  by  Ai'(-x)  and  Bi'(-x),  respectively. 

H.  MODE  DEPTH  FUNCTIONS 

The  normalized  mode  depth  functions  are  given  by 

Un(Z0)  Un(Z)  =  Fn(Z0  Fn(Z)/Dn  (82) 

where  n  is  the  mode  number,  Z0  is  the  source  depth,  and  Z  is  the  receiver  depth.  The 
Fn(Z0)  or  Fn(Z)  is  the  Fj  of  eq.  72  for  the  layer  in  which  the  source  or  receiver 
respectively  are  located,  as  evaluated  for  mode  n. 

The  major  objective  of  this  section  is  to  provide  expressions  for  Dn,  the  normaliza¬ 
tion  factor.  Reference  1  defines  a  canonical  normalization  factor  Dc  related  to  Dn  by 

Dn  =  af 1  Dc  .  (83) 

The  significant  feature  of  Dc  is  that  it  can  be  expressed  in  terms  of  canonical  vari¬ 
ables. 

However,  in  the  preparation  of  this  report,  we  found  an  improved  method  of  treat¬ 
ment.  If  we  substitute  eq.  53  for  the  frequency  dependence  in  a-|,  eq.  83  may  be  writ¬ 
ten  as 

Dn  =  J  Dca  (84) 
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where 


Dca  =  |y  -  x|_1|Dc| 


(85) 


and 


J  =  Ci|Ci2  C2~2  -1|  I2/J-1  .  (86) 

Equation  85  defines  a  new  canonical  normalization  factor  that  can  be  converted  to  Dn 
by  the  scale  factor  of  eq.  86.  The  original  normalization  factor  required  two  plots  such 
as  presented  in  reference  3.  The  first  was  a  plot  of  Dc  versus  x  and  the  second  was  Dn 
versus  frequency.  However,  these  two  plots  may  be  replaced  by  a  single  plot  in  which 
we  plot  Dca  versus  fca-  The  application  of  the  scale  factor  J  to  the  ordinate  and  the 
scale  factor  K  to  the  abscissa  gives  Dn  as  a  function  of  frequency.  Numerical  examples 
of  this  will  be  presented  later. 

If  one  solves  eq.  1  for  Z2,  one  obtains  the  right  side  of  eq.  86.  Thus,  the  scale 
factor  J  is  the  thickness  of  layer  1.  This  was  not  recognized  in  the  original  deviation 
until  we  observed  that  the  numerical  values  of  J  and  Z2  were  the  same  for  a  particular 
example. 

The  next  step  is  to  present  the  various  expressions  for  Dc  for  use  in  eq.  85.  From 
eq.  90  of  reference  1,  Dc  may  be  expressed  as 


Dc  =  Ds  +  Db  +  D„i  (87) 

where  Ds  and  Db  are  contributions  at  the  upper  and  lower  boundaries  respectively  and 
Dn  is  a  contribution  at  the  interface  between  layers  1  and  2.  From  eqs.  93  and  97  of 
reference  1 

Dni  =  (1 +p3)[y  Fj2(-y)+  Fi'2(-y)]  .  (88) 

Here,  F,  and  F/  are  given  by  eqs.  77  and  79  or  their  rigid  surface  counterparts. 

In  eq.  87,  Ds  and  Db  depend  on  the  surface  and  bottom  boundary  conditions.  For  a 
free  or  rigid  surface, 

Ds  =  -tt  "2  (89) 

or 

Ds  =  -  xtr-2  .  (90) 
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For  a  free  or  rigid  bottom, 
Db  =  -P  F2'2(-£20) 


(91) 


or 


Db  =  -  P  £20  F22  (~  V2o) 

(92) 

where 

F2'(-  C20)  =  -  D2  Ai'(-  £20)  -  E2Bi'(-  £20) 

(93) 

F2(-  £20)  =  D2  Ai(-  £ 20 )  +E2Bi'(-  £20) 

(94) 

and  £20  is  given  by  eq.  11.  Equations  89  to  94  follow  respectively  from  eqs.  94,  95, 

98,  99,  76,  and  73  of  reference  1. 

I.  DETERMINATION  OF  EIGENVALUES  FOR  A  GIVEN  FREQUENCY 

Consider  the  case  where  one  wants  to  determine  the  eigenvalues  for  a  given  fre¬ 
quency,  f0.  For  example,  we  may  wish  to  sum  up  the  modes  at  a  given  frequency  to 
determine  the  field  or  we  may  wish  to  compare  the  standing  wave  patterns  of  modes  at 
a  given  frequency.  This  requires  a  special  computer  routine  because  the  points  at 
which  the  canonical  eigenvalues  are  available  will  lead  to  different  frequencies  for  the 
different  modes.  We  cannot  rely  upon  some  simple  interpolation  routine  because  the 
sensitivity  of  the  mode  evaluation  at  or  near  the  critical  frequencies  require  a  high 
degree  of  accuracy. 

From  eqs.  53  and  54,  we  note  that 

f„=|yo-x„|3/2K.  (95) 

Here,  f0  is  the  given  frequency  and  K  is  known  from  the  profile  parameters.  Our  task 
is  to  determine  y0  and  x0.  From  eq.  95  it  follows  that 

y0  =  x0  +  (Kfo1)273  (96) 

We  now  solve  the  eigenvalue  equation 


G(x,  y)  =  0  (97) 

under  the  constraint  of  eq.  96.  We  iterate  eq.  97  to  determine  x0  by  substituting  eq.  96 
for  y  in  eq.  97.  The  derivative  of  eq.  97  is  given  by 

dG/dx  =  <5G/<5x  +  (<3G/<5y)  (dy0/dxo)  .  (98) 
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From  eq.  96  it  follows  that  dy0/dx0  =  1.  Thus,  eq.  98  reduces  to 

dG/dx  =  <5G/<5x  +  <5G/<5y  .  (99) 

Starting  estimates  of  are  obtained  from  the  plots  of  canonical  frequency  ver¬ 
sus  x.  The  expressions  for  G(x,y),  8G/8x,  and  SG/8y,  appropriate  to  the  boundary 
conditions,  are  used  in  the  iteration,  where  eq.  96  is  always  substituted  for  y  through¬ 
out  the  evaluation. 

J.  ESTIMATES  OF  EIGENVALUE  DIFFERENCES  AT  THE  CRITICAL 
FREQUENCIES 

Section  2F  presented  estimates  of  the  critical  frequencies  in  terms  of  the  parame¬ 
ters,  S,  p,  and  Airy  function  zeros.  In  this  section,  we  present  approximations  to  the 
eigenvalue  differences  at  the  estimate  of  critical  frequencies.  These  differences  are 
important  because  they  are  a  measure  of  the  amount  of  coupling  between  the  two 
modes  associated  with  the  critical  frequency. 

Appendix  A  derives  expression  for  approximate  eigenvalues  based  on  a  Taylor 
series  expression  about  x0  for  fixed  y0.  Here.  \u  und  y0  represent  a  point  of  intersec¬ 
tion  of  the  vertical  asymptotes  of  ens.  4C  or  41  with  the  slant  asymptotes  of  eqs.  42 
or  44.  This  expression  is 

(x-x0)2  +  b(x-x0) +  c  =  0  .  (100) 

The  two  values  of  x  for  this  quadratic  equation  represent  an  approximation  to  the 
eigenvalues  for  the  two  modes  associated  with  the  intersection  of  asymptotes.  These 
values  of  x  correspond  to  y  =  y0. 

First,  consider  the  case  of  the  free-rigid  boundary  condition.  Here, 
b  =  KnP-2(l  -S)-,(G2/G4)  +  MG3/G4)  (101) 

and 

c  =  Mmnp-2(l-S)-,(G,/G4)  .  (102) 

The  coefficients  Kn,  L,,,,  and  Mmn  are  constants  that  depend  upon  the  particular  inter¬ 
section  point.  The  coefficients  have  no  connection  with  the  factors  K  and  L  previously 
given.  They  are  given  by 


K„  =  -  [Tra'n  Ai2(a'n)]_1 
Lm  =  -  [7rAi'2(am)]_1 
Mmn  =  Ljn  Kn  . 


(103) 

(104) 

(105) 


Here,  am  is  the  negative  of  the  m'th  zero  of  Ai  and  a'n  is  the  negative  of  the  n'th  zero 
of  Ai'.  Here, 


X0  — 


(106) 


and 


y0  =  (1  -  S_1)am  +  a'n/p2S  .  (107) 

For  the  free-free  boundary  condition,  eqs.  101,  102,  104,  and  105  apply  with 
Kn  =  Ln  (108) 

and 

Mmn  =  LmLn.  (109) 

Equation  107  is  replaced  by 

y0  =  (1  —  S_1)am  +  an/p2S  .  (110) 

For  the  rigid-rigid  boundary  condition,  eqs.  101,  102,  and  103  apply  with 


L^K 


m 


(in) 


and 


Mmn  =  KmKn.  (112) 

Equations  106  and  107  are  replaced  by 

x0  =  a'n  .  (113) 


and 


yQ  =  (1  -S"1)a'm  +  a'n/p2S  . 


(114) 
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K.  PHASE  VELOCITY  DIFFERENCES  AND  MODE  INTERFERENCE 
RANGES 

First,  we  develop  an  expression  for  the  difference  in  phase  velocity  between  the 
adjacent  modes  at  the  critical  frequency.  Let  xn+],  yn+i  and  xn,  yn  be  the  eigenvalues 
for  mode  n  +  1  and  mode  n  at  the  critical  frequency.  If  one  expands  eq.  55  by  the 
binomial  expression  for  mode  n  +  1  and  for  mode  n  and  subtracts,  one  obtains 

CPi„+i-Cp>n  =  C1r2/3ir-2/3|y1|2/32'1(xn+i-xn)  .  (115) 

The  substitution  of  eq.  53  into  eq.  115  leads  to 

Cp.n+1  -  Cp,n  =  C2(l  -  Cj2  C2'2  )2'1(x  -  y)',(xn+1  -  xn)  .  (116) 

Since  the  frequency  is  the  same  for  both  modes,  x  -  y  may  be  evaluated  from  either 
eigenvalue. 

Estimates  of  the  phase  velocity  differences  corresponding  to  the  estimates  of  sec¬ 
tion  2J  are  readily  obtained  with  the  use  of  eq.  116.  Here,  (xn+1  -xn)  is  given  by  the 
discriminate  of  eq.  100.  In  the  case  of  the  free-rigid  boundary  condition,  x  and  y  are 
estimated  by  eqs.  106  and  107.  In  the  case  of  the  free-free  boundary  condition,  y  is 
estimated  from  eq.  107;  and  in  the  case  of  the  rigid-rigid  boundary  condition,  x  and  y 
are  estimated  from  eqs.  113  and  114. 

Reference  4  discusses  the  mode  interference  distance.  Since  modes  have  slightly 
different  phase  velocities,  the  phase  relationship  between  modes  changes  with  range. 
The  mode  interference  distances  is  the  difference  between  successive  range  points 
where  the  modes  are  in  (or  out)  of  phase.  The  mode  interference  distance  may  be 
expressed  as 

Rn+l  =  f  1  Cp  n  Cp  nti  [Cp,n+i  -  Cpn]  1  .  (f  17) 

If  we  substitute  eqs.  53  and  116  into  eq.  117  and  simplify,  we  obtain 

Rn+i  =  2C2'ir|y1|~i|Ci2  C2~2  -  l|1/2(x -y)'1/2(xn+1 -x„)_1  .  (118) 

One  may  estimate  Rn+i,  based  on  the  estimates  of  section  2J,  as  just  discussed  for  esti¬ 
mates  of  the  phase  velocity  differences. 
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3.  NUMERICAL  EXAMPLES 


This  section  presents  numerical  examples  of  two  different  types  of  functions,  i.e., 
the  canonical  and  physical  functions.  The  canonical  functions  are  those  which  are 
completely  determined  by  the  parameters  p  and  S  and  the  boundary  conditions.  These 
include  the  canonical  eigenvalues,  the  canonical  frequency,  the  canonical  group  velocity 
factor,  and  canonical  eigenfunction  normalization  coefficient.  The  physical  functions 
are  those  which  are  determined  from  the  canonical  functions  and  a  particular  choice  of 
the  sound-speed  profile  parameters.  The  first  of  these  functions  relates  acoustic  fre¬ 
quency  to  the  eigenvalue  solutions.  Three  other  physical  functions  are  plots  of  phase 
velocity,  group  velocity,  and  eigenfunction  normalization  coefficients,  all  given  as  func¬ 
tions  of  frequency.  The  last  of  the  physical  functions  are  plots  of  the  eigenfunctions. 
Here,  the  eigenfunctions  (standing  waves)  are  plotted  as  a  function  of  receiver  depth 
for  a  particular  choice  of  frequency  and  source  depth. 

A.  FREE-RIGID  BOUNDARY,  p  =  1,  S  =  0 

This  particular  case  was  chosen  for  first  presentation  because  this  case  is  the  sim¬ 
plest  double  duct  with  no  coupling  effects  between  the  ducts.  For  S  =  0,  the  surface 
and  bottom  sound  speeds  of  figure  1  are  the  same.  For  p  =  1,  the  slopes  at  the  inter¬ 
face  between  layers  are  equal  but  of  opposite  sign.  Thus,  this  is  the  class  of  all  two- 
layer  profiles  that  are  symmetric  about  an  axis  of  minimum  or  maximum  sound  speed 
and  with  free-rigid  boundary  conditions. 

Figure  2  presents  the  canonical  eigenvalues  for  this  case.  We  point  out  some  fea¬ 
tures  common  to  all  canonical  eigenvalue  plots  of  this  report.  Each  plot  contains  one 
set  of  horizontal  asymptotes.  These  are  the  solid  horizontal  grid  lines,  and  are  solu¬ 
tions  to  eq.  46.  For  the  case  of  p  =  1,  this  reduces  to  y  =  an  or  a'n.  Each  plot  contains 
two  sets  of  asymptotes  for  y  =  -  »,  shown  by  the  vertical  dotted  lines.  One  set  of 
these  asymptotes  is  always  vertical  and  independent  of  the  values  of  p  and  S.  For  the 
free  surface,  these  asymptotes  are  given  by  x  =  an .  The  second  set  of  asymptotes  is 
usually  a  slant  set  except  for  the  special  case  of  S  =  0.  Here,  they  are  vertical  and  for 
p  =  1,  the  solution  for  a  rigid  bottom  is  x  =  a'n.  Since  the  critical  frequencies  are 
associated  with  the  intersection  of  the  vertical  and  slant  asymptotes,  there  are  no  criti¬ 
cal  frequencies  for  S  =  0  because  the  slant  asymptote  became  vertical  for  this  special 
case. 

Next,  we  point  out  various  properties  that  apply  to  all  of  the  eigenvalue  plots  of 
this  report.  These  properties  also  apply  to  figures  2  to  4  of  reference  3  for  single 
ducts.  Indeed,  these  properties  are  fundamental  to  the  canonical  approach.  We  state 
these  properties  without  justification  because  such  justification  may  be  found  in  refer¬ 
ence  3. 
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Figure  2.  Eigenvalues  for  S  =  0,  p  =  1,  free-rigid  boundary. 

Observe  the  line  y  =  x.  This  line  is  given  in  all  eigenvalue  plots  because  it  divides 
the  modes  into  two  families,  with  the  modes  below  and  above  the  line  associated  with 
cases  A  and  B  respectively  of  figure  1.  The  mode  count  increases  as  each  family  mem¬ 
ber  moves  away  from  the  line  y  =  x.  To  be  more  specific,  the  modes  are  ordered  by 
their  displacement  from  the  line  y  =  x  for  fixed  y.  Unless  otherwise  specified,  we  will 
present  the  first  five  modes  of  each  family.  The  same  method  of  coding  will  be  used 
to  associate  the  eigenvalues  and  other  curves  with  a  particular  mode  number.  Mode  1 
curves  will  be  solid.  Mode  2  curves  will  consist  of  long  dashes  separated  by  a  dot. 
Mode  3  curves  will  consist  of  short  dashes,  Mode  4  curves  will  consist  of  short  and 
long  dashes.  Mode  5  curves  will  consist  of  long  dashes  separated  by  two  dots. 

The  large  negative  values  of  y  or  x  correspond  to  high  frequencies.  The  frequency 
decreases  as  we  move  along  a  given  case  A  mode  to  larger  values  of  y  or  along  a 
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given  case  B  mode  to  larger  values  of  x.  For  large  positive  values  of  y  or  x,  the  curves 
approach  the  line  y  =  x  as  an  asymptote.  These  values  correspond  to  low  frequencies 
because  the  frequency  is  zero  for  the  line  y  =  x. 

For  the  A  set  of  modes,  the  point  where  the  curves  cross  y  =  0  corresponds  to  the 
ray  with  phase  velocity  equal  to  the  interface  sound  speed.  Thus,  the  eigenvalues  with 
positive  y  correspond  to  rays  reflecting  from  both  boundaries.  The  eigenvalues  with 
negative  y  correspond  to  modes  trapped  in  the  upper  or  lower  duct.  In  this  case,  the 
odd  modes  correspond  to  rays  reflecting  from  the  bottom  and  form  a  refractive  apex 
in  layer  2  while  the  even  modes  correspond  to  rays  reflecting  from  the  surface  and 
form  a  refractive  nadir  in  layer  1.  This  association  is  made  because  the  odd  asymp¬ 
totes  are  determined  by  x  =  a'n  for  the  rigid  bottom  while  the  even  asymptotes  are 
determined  by  x  =  a„  for  the  free  surface. 

For  the  B  set  of  modes,  the  point  where  the  curves  cross  x  =  0  corresponds  to  the 
ray  with  phase  velocity  equal  to  the  surface  sound  speed.  Thus,  positive  values  of  x 
correspond  to  modes  reflecting  from  the  ocean  surface  (and  also  from  the  bottom  for 
this  case  of  S  =  0),  while  negative  values  of  x  correspond  to  rays  that  form  an  apex  in 
layer  1.  For  this  case  of  S  =  0,  they  also  form  a  nadir  in  layer  2. 

Note  that  the  spacing  between  the  eigenvalues  of  a  given  family  is  relatively  uni¬ 
form.  This  is  not  surprising  for  the  case  B  modes  because  this  is  a  characteristic  prop¬ 
erty  of  all  single  ducts.  As  we  shall  see  presently,  this  is  also  true  for  the  case  A 
modes  reflecting  from  both  boundaries  (y  >  0),  but  is  not  generally  true  for  case  A 
modes  for  y  <  0.  We  note  that  for  y  somewhat  greater  then  0,  the  spacing  between 
mode  1  of  case  A  and  mode  1  of  case  B  is  also  comparable  to  the  spacing  between 
modes  1  and  2  of  cases  A  or  B.  This  property  will  take  on  more  significance  when  we 
investigate  the  other  two-boundary  conditions. 

The  case  of  figure  2  is  unique  in  that  the  cases  A  and  B  eigenvalues  are  symmetric 
about  the  line  y  =  x.  This  may  be  demonstrated  by  simplifying  eq.  28  for  the  case  of 
p  =  1  and  S  =  0,  and  noting  that  it  remains  the  same  eigenvalue  equation  when  x  or  y 
are  interchanged. 

The  same  general  approach  was  used  to  generate  all  eigenvalue  curves  of  this 
report,  such  as  those  of  figure  2.  Newton’s  method  was  applied  using  eq.  29  (or  eq.  32 
or  35  for  the  other  boundary  conditions)  for  fixed  values  of  y  for  the  case  A  modes. 
Equation  30  (also  33  or  36)  was  used  for  fixed  values  of  y  for  the  case  B  modes.  The 
reason  for  this  choice  is  that  Newton’s  method,  using  eq.  29,  fails  for  the  case  A 
modes  when  the  slope  is  vertical  and  Newton’s  method,  using  eq.  30,  fails  for  the 
case  B  modes  when  the  slope  is  horizontal.  The  iterative  process  for  the  cases  A  and  B 
modes  is  commenced  at  large  negative  values  of  y  and  x  respectively  as  given  by  the 
asymptotic  solutions. 


25 


Although  we  will  present  three  other  canonical  functions  that  depend  only  on  p  and 
S  and  are  independent  of  the  specific  profile  parameters,  it  is  onvenient  to  introduce 
specific  parameters  at  this  point.  For  all  of  the  case  A  examples  of  this  report,  we 
choose  C,  =  1480  m/s,  C2  =  1482  m/s,  and  7i  =  0.02  s'1.  These  are  the  same  parame¬ 
ters  used  in  reference  3  for  one-layer  ducts.  From  eq.  1,  we  determine  that 
Z2  =  99.797661  m  and  from  eq.  3  that  7i0  =  0.020081186  s*1.  This  gives  all  of  the 
useful  parameters  of  the  first  layer.  The  parameters  of  the  second  layer  may  then  be 
determined  from  these  first  layer  parameters  and  the  values  of  p  and  S. 

For  our  values  of  p  and  S,  it  follows  that  C3  =  C1t  Z3  =  2  Zs,  72  =  -7]0,  and  720  = 
-7].  The  left  panel  of  figure  3  is  a  plot  of  sound-speed  profile  used  for  the  case  A 
modes. 

The  left  scale  ot  figure  4  presents  the  canonical  frequency,  plotted  as  a  function  of 
x.  The  right  scale  gives  the  acoustic  frequency  as  determined  by  eq.  84.  The  scale 
factor  J  is  given  by  eq.  86  and  is  equal  to  45.45  for  the  profile  of  figure  3.  This  same 
value  applies  to  all  of  the  case  A  modes  of  this  report. 

Figure  5  presents  the  phase  velocity  as  determined  by  eq.  56.  Figure  5  displays  a 
typical  behavior  for  a  single  duct  with  a  relatively  uniform  spacing  between  modes.  At 
high  frequencies,  all  curves  asymptotically  approach  1480  m/s,  the  sound  speed  at  the 
two  boundaries.  At  low  frequencies,  the  curves  of  figure  5  form  vertical  asymptotes  at 
the  cut-off  frequencies  of  the  duct.  The  evaluation  of  these  cut-off  frequencies  is  very 
involved  and  will  not  be  attempted  here.  Reference  3  presents  such  an  evaluation  for  a 
one-layer  duct. 

Figure  6  presents  the  canonical  group  velocity  factor  of  eq.  58  as  evaluated  using 
eq.  59.  Figure  7  presents  the  group  velocity  as  a  function  of  frequency.  This  was 
evaluated  by  eq.  57  with  the  use  of  L  of  figure  6.  At  high  frequencies,  all  curves 
asymptotically  approach  1480  m/s,  the  sound  speed  at  the  two  boundaries.  As  the 
frequency  is  decreased  to  low  values,  all  curves  go  rapidly  toward  small  values  of 
group  velocity  corresponding  to  steep  rays  reflecting  off  both  boundaries.  At  intermedi¬ 
ate  frequencies,  the  group  velocity  forms  a  relative  maximum  (point  of  Airy  phase). 
Our  previous  work  with  a  positive-gradient  surface  duct  overlaying  a  negative-gradient 
half  space  yields  results  very  similar  to  this.  We  found  that  at  high  frequencies  the 
relative  maximum  goes  to  the  ray  theory  value  for  the  ray  which  grazes  the  bottom  of 
the  surface  layer.  This  is  a  group  velocity  one-third  of  the  distance  between  the  surface 
sound  speed  and  the  sound  speed  at  the  bottom  of  the  surface  layer.  For  the  profile  of 
figure  3,  this  is  a  value  of  1480.67  m/s,  plotted  as  a  horizontal  line  in  figure  7.  We  see 
that  the  same  limit  appears  to  apply  here. 
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Figure  5.  Phase  velocities  for  the  A  modes 
of  figure  3. 
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Figure  6.  Group  velocity  factor  for  the  A  modes  of  figure  2 
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Figure  7.  Group  velocities  for  the  A 
modes  of  figure  3. 


The  left  scale  of  figure  8  gives  the  canonical  normalization  factor  of  eq.  85  as 
evaluated  by  eq.  87  to  94.  The  right  scale  gives  the  customary  eigenvalue  normaliza¬ 
tion  coefficient  as  defined  by  eq.  84.  The  factor  J  is  the  thickness  of  layer  1  or  99.8  m 
for  the  profile  of  figure  3.  The  bottom  scale  is  the  canonical  frequency,  whereas  the 
top  scale  is  the  acoustic  frequency  of  eq.  53.  The  curves  for  modes  2  and  4  (modes 
trapped  in  the  upper  duct)  are  well  behaved.  However,  the  curves  for  modes  1,  3, 
and  5  rapidly  grow  rather  large.  The  dominant  terms  of  Dc  lie  in  eq.  88  and  are  y 
Bi2(-y)  and  Bi'2(-y).  Both  terms  diverge  as  TT'^-y)172  exp  [4(-y)3/2/3]  for  large  nega¬ 
tive  y,  i.e.,  for  large  frequency.  The  coefficient  for  these  terms  is  Ai2(-x)  as  can  be 
seen  from  eqs.  77  and  79.  This  coefficient  is  zero  for  the  even  modes  and  the  diver¬ 
gent  terms  drop  out.  However,  the  coefficient  is  not  zero  for  the  even  modes  and  Dc 
diverges  for  high  frequency.  We  note  that  for  the  case  of  p  =  1,  eqs.  80  and  81  reduce 
to  D2  =  Bi(-x)  and  Ej  =  Ai(-x).  Thus,  the  term  of  eq.  92  poses  no  problem. 

As  we  shall  see  later,  there  are  cases  for  which  the  Dn  for  odd  modes  do  not 
diverge,  but  form  a  relative  maximum.  From  eq.  85  we  note  that  even  though  Dc 
diverges,  Dca  and  Dn  could  converge  because  (y-x)'1  goes  to  zero  in  the  high-frequency 
limit.  We  evaluated  the  Dca  of  figure  8  using  double-precision  routines.  However,  the 
Dca  for  odd  modes  appeared  to  diverge  rather  than  form  a  relative  maximum. 
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Figure  8.  Canonical  or  customary 
eigenvalue  normalization  coefficients 
for  the  A  modes  of  figures  2  or  3. 

The  physical  explanation  for  the  different  behavior  of  Dn  for  even  and  odd  modes 
lies  in  the  fact  that  in  eq.  76  we  have  chosen  the  unnormalized  depth  function  to  sat¬ 
isfy  the  boundary  condition  at  the  surface.  This  is  not  a  good  choice  for  modes  highly 
trapped  in  the  lower  duct.  We  could  choose  the  unnormalized  depth  function  to  satisfy 
the  bottom  boundary  condition.  This  was  the  technique  used  to  resolve  a  similar  prob¬ 
lem  for  the  single-layer  ducts  of  reference  3.  However,  in  this  approach  the  even 
modes  would  then  diverge. 

Huge  or  divergent  values  of  Dn  do  not  pose  a  problem  mathematically,  because  Dn 
grows  large  with  increasing  frequency  to  compensate  for  the  large  growth  in  the  unnor¬ 
malized  depth  functions.  The  only  problem  this  poses  is  in  the  presentation  of  Dn  for 
both  upper  and  lower  duct  modes  on  the  same  scale. 

We  are  not  completely  satisfied  with  our  present  approach,  which  yields  too  large  a 
range  of  values  for  Dn.  However,  one  can  multiply  the  unnormalized  eigenfunction  by 
any  constant  (i.e.,  a  value  independent  of  depth)  and  this  constant  will  appear  squared 
in  Dn  to  yield  the  normalized  depth  function.  Thus,  Dn  is  by  no  means  unique  but 
depends  on  how  the  unnormalized  eigenfunctions  are  defined.  For  example,  an 
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alternative  approach  to  that  illustrated  in  figure  8  is  to  divide  the  present  unnormalized 
eigenfunction  of  eq.  72  for  odd  modes  by  the  term  Bi'  (-y)  and  carry  out  the  normali¬ 
zation  process  in  the  same  manner  as  before. 

We  will  point  out  some  other  interesting  properties  of  the  normalization  coefficient 
after  we  have  presented  other  numerical  examples.  There  appear  to  be  many  aspects 
of  Dn  that  are  of  further  interest.  However,  the  scope  of  the  present  report  will  not 
permit  a  thorough  analysis. 

Figure  9  presents  the  normalized  depth  functions  at  500  Hz  for  A  modes  1  to  5  as 
a  function  of  receiver  depth  for  a  fixed  source  depth  for  a  given  mode.  These  curves 
are  generated  from  eq.  82  with  the  use  of  the  appropriate  normalization  values,  associ¬ 
ated  with  figure  8.  The  eigenvalues  at  500  Hz,  used  in  this  evaluation,  were  deter¬ 
mined  by  the  approach  discussed  in  section  21.  The  source  depth  is  170,  27,  153,  60, 
and  125  m  for  modes  1  to  5  respectively.  These  depths  were  chosen  to  enhance  the 
particular  mode  so  that  the  behavior  of  the  curve  for  small  amplitudes  could  be  deter¬ 
mined.  The  dotted  horizontal  line  represents  the  interface  between  layers  at  99.8  m. 
The  bottom  of  the  channel  is  at  twice  this  depth  or  199.6  m.  The  dotted  vertical  line  is 
the  zero  origin  for  the  depth  functions. 

All  of  the  curves  go  to  zero  at  the  surface  and  have  zero  slope  (form  an  end  point 
maximum)  at  the  bottom. 

There  is  one  method  of  mode  identification  from  depth  function  plots  that  always 
gives  the  correct  mode  number.  The  mode  number  is  one  more  than  the  number  of 
zero  crossings.  This  rule  obviously  gives  the  mode  number  for  modes  3  to  5,  but  it  is 
not  quite  so  obvious  for  modes  1  and  2.  However,  an  examination  of  detailed  com¬ 
puter  printouts  reveals  that  mode  1  values  are  always  positive  and  that  mode  2  values 
turn  negative  toward  the  bottom  of  layer  2.  In  this  particular  example,  the  mode  num¬ 
ber  is  also  equal  to  the  number  of  nodes  (zeros)  or  the  number  of  antinodes  (relative 
or  end  point  maxima).  However,  later  we  will  present  other  examples  where  these 
criteria  for  mode  identification  fail. 

Recall  that  in  the  discussion  of  figure  2,  we  associated  odd  modes  with  the  lower 
duct  and  even  modes  with  the  upper  duct.  In  figure  9,  modes  1  and  3  are  clearly 
trapped  in  the  lower  duct  and  modes  2  and  4  in  the  upper  duct.  Mode  5  has  a  phase 
velocity  slightly  above  1482  m  and  is  not  trapped  in  either  duct,  although  the  ampli¬ 
tudes  are  somewhat  larger  in  the  lower  layer.  Later  on,  we  will  see  some  more  impres¬ 
sive  examples  of  how  the  eigenvalue  plots  may  be  used  to  predict  the  duct  in  which 
the  mode  is  trapped,  i.e.,  the  mode  has  large  eigenfunction  amplitudes. 
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Figure  9.  Eigenfunctions  at  500  Hz  for  the  A  modes  of  figure  3. 
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For  the  most  part,  we  have  little  interest  in  the  case  B  modes  because  they  repre¬ 
sent  propagation  in  a  bounded  refractive  duct  with  none  of  the  characteristics  of  dou¬ 
ble  duct  propagation.  However,  the  symmetry  about  the  line  y  =  x  in  figure  2  makes 
the  case  B  modes  of  interest  for  this  special  case.  We  evaluate  these  modes  for  the 
profile  parameters  given  by  C,  =  1482  m/s,  C2  =  1480  m/s,  and  =  -0.02(C,/C2)3  = 
-0.020081186  s'1.  This  is  the  profile  shown  in  the  right  panel  of  figure  3.  These 
parameters  are  chosen  so  as  to  have  a  particular  relationship  with  the  parameters  of 
the  left  panel.  Let  the  unbarred  values  be  the  parameters  for  the  case  A  profile  and 
barred  values  be  those  for  the  case  B  profile.  Then,  the  following  relationship  holds: 

Ci  =  C2,  C2  =  Cj,  C3  =  C2,  yx  =  - yi0,  y]0  =  - y,,  ¥2  =  Y}  and  y2Q  =  y]Q.  Furthermore,  from 
the  symmetry  properties  x  =  y  and  y  =  x.  It  follows  then  from  eq.  53  to  55  that 

fca  =  fca,  K  =  K,  and  f  =  f  , 

i.e.,  the  frequency  corresponding  to  the  point  x,  y  for  case  A  is  the  same  as  that  to  the 
point  x,  y  for  case  B. 

Figure  10  presents  the  counterpart  of  figure  4  for  the  case  B  modes.  The  dependent 
variable  is  the  same  in  figure  10  as  in  figure  4.  Indeed,  if  we  had  plotted  the  depend¬ 
ent  variable  as  a  function  of  y  rather  than  x  in  figure  10,  it  would  be  identical  to  the 
plots  of  figur:  4. 


Figure  10.  Canonical  or  acoustic  frequencies 
for  the  B  modes  of  figures  2  or  3. 
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The  phase  velocity  counterpart  of  figure  5  for  the  case  B  modes  was  evaluated  and 
is  identical  to  figure  5.  To  demonstrate  this,  we  compare  the  result  of  eq.  56  as  evalu¬ 
ated  for  the  B  modes  with  that  of  eq.  30  of  reference  3  as  evaluated  for  the  A  modes. 
These  expressions  are  identical  because  Q  =  C2,  f  =  f,  l/J  =  |y10l>  and  x  =  y  . 

Figure  11  is  the  counterpart  of  figure  6  for  the  case  B  modes.  The  results  of  fig¬ 
ure  11  were  used  to  evaluate  the  group  velocity  for  the  case  B  modes.  This  plot  was 
identical  to  that  of  figure  7  for  the  A  modes.  This  property  could  be  demonstrated  by 
comparing  the  two  evaluations  of  eq.  57.  However  a  much  easier  demonstration  fol¬ 
lows  from  eq.  63  of  reference  1.  Since  the  phase  velocity  is  the  identical  function  of 
frequency  for  the  two  cases,  eq.  63  easily  demonstrates  that  the  group  velocities  are 
also  identical  functions  of  frequency. 


x 

Figure  11.  Group  velocity  factor  for  the  B  modes  of  figure  3. 


Figure  12  is  the  B  mode  counterpart  of  figure  8.  In  contrast  to  figure  8,  Dn  diverges 
with  increasing  frequency  for  both  even  and  odd  modes.  Here  the  dominant  terms  of 
Dn,  in  eq.  88  involves  Bi'2(-x).  The  coefficient  of  this  term  is  Ai'2(-y)  in  eq.  77  and 
[Ai'(-y)j2  in  eq.  79.  Here,  Ai(-y)  and  Ai'(-y)  cannot  both  be  zero.  Thus,  Dc  diverges 
for  both  A  and  B  modes  for  large  frequency.  We  note  that  Dc  of  eq.  92  reduces  to  a 
term  containing  -x  Bi2(-x)  Ai?(-x),  that  has  a  finite  limit  for  large  negative  x  and 
poses  no  problem. 
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Figure  12.  Canonical  or  customary  eigenvalue 
normalization  coefficients  for  the  B  modes  of 
figures  2  or  3. 

Two  other  approaches  might  result  in  a  nondivergent  form  for  Dn.  One  is  to  choose 
the  unnormalized  depth  functions  so  as  to  automatically  satisfy  one  of  the  continuity 
conditions  at  the  axis.  The  other  approach  is  to  divide  the  eigenfunctions  of  eq.  72  by 
the  term  Bj(-x).  A  similar  approach  has  already  been  discussed  for  the  A  modes. 

Figure  13  is  the  B  counterpart  of  figure  9.  The  source  depths  are  100,  75,  60,  40, 
and  30  m  for  modes  1  to  5  respectively.  In  contrast  to  case  A,  the  mode  amplitudes 
are  comparable  in  the  upper  and  lower  layers.  This  is  the  case  of  a  single  refractive 
duct  with  symmetric  sound  speeds  about  the  axis.  The  odd  modes  have  a  node  near 
the  axis  while  the  even  modes  have  an  antinode  near  the  axis.  These  features  do  not 
occur  exactly  at  the  axis  because  the  free-rigid  boundary  conditions  break  the  symme¬ 
try. 

To  further  compare  the  results  of  modes  A  and  B,  we  note  from  eq.  75  that  al  for 
the  B  modes  is  -aj  for  the  A  modes.  Moreover,  it  follows  from  eq.  73  that  =  £2. 
This  means  that  the  Airy  function  arguments  in  layer  2  of  figure  9,  as  a  function  of 
depth  from  the  axis,  is  the  same  as  these  arguments  in  layer  1  of  figure  13  as  a  func¬ 
tion  of  depth  from  the  surface.  This  property  shows  up  for  mode  1  in  figures  9  and  13 
but  to  little  degree  in  modes  2  to  5.  There  are  two  reasons  for  this  lack  of  comparison. 
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Figure  13.  Eigenfunctions  at  500  Hz  for  the  B  modes  of  figure  3. 
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The  most  important  reason  is  that  the  coefficients  Ds  and  Ej  are  different  due  to  differ¬ 
ent  requirements  at  the  surface  and  axis.  These  coefficients  weight  the  contribution  of 
Ai(-£)  and  Bi(-£)  in  a  different  manner  and  produce  a  different  response  even  though 
£  is  the  same.  A  second  reason  is  that  the  source  depths  are  not  the  same  for  the  A 
and  B  modes.  This  scales  the  amplitudes  by  a  multiplication  factor,  but  would  not 
affect  the  shape  as  a  function  of  depth. 

One  may  also  show  that  ^  Here,  even  mode  1  in  the  upper  layer  of  figure  9 
shows  little  resemblance  to  mode  1  in  the  lower  layer  of  figure  13  because  of  radically 
different  weighting  of  D;  and  Ej. 

The  important  feature  to  note  is  that  the  eigenfunctions  are  radically  different  in 
figure  9  than  those  of  figure  13.  This  occurs  even  though  the  eigenvalues  are  identical. 

Our  interest  in  comparing  the  cases  A  and  B  results  stems  from  a  classical  mathe¬ 
matics  problem  that  was  the  subject  of  a  film  produced  by  the  American  Mathematical 
Society  entitled  “Can  You  Hear  the  Shape  of  a  Drum?”  The  film  discussed  the  extent 
to  which  the  environment  could  be  inferred  from  the  eigenvalues.  This  is  related  to  the 
inversion  problem  in  which  physical  properties  of  the  environment  are  determined  from 
acoustic  field  measurements  of  travel  time. 

Our  example  demonstrates  that  the  eigenvalues  (phase  velocity  versus  frequency) 
do  not  uniquely  determine  the  environment  nor  do  the  group  velocity  characteristics. 

On  the  other  hand,  the  eigenfunctions  (pressure  field  versus  depth)  are  markedly  dif¬ 
ferent  for  the  two  environments.  The  problem  of  whether  the  eigenfunctions  uniquely 
determine  the  environment  will  not  be  addressed  here. 

Our  example  of  two  radically  different  environments  that  produce  identical  group 
velocity  characteristics  suggests  that  measurements  of  travel  time  differences  alone  may 
not  suffice  to  uniquely  define  the  environment.  This  would  be  most  likely  in  the  case 
of  long-range  measurements  where  travel  time  differences  are  dominated  by  the  group 
velocity  characteristics. 

B.  FREE-RIGID  BOUNDARY,  p  =  1.2,  S  =  0.5 

This  case  was  chosen  to  illustrate  the  coupling  phenomenon  between  double  ducts. 
This  case  was  presented  at  the  120th  Meeting  of  the  Acoustical  Society  of  America 
described  in  reference  7. 

Figure  14  presents  the  canonical  eigenvalues.  The  B  modes  have  horizontal  asymp¬ 
totes  given  by  the  solution  to  eq.  46.  The  A  modes  have  two  sets  of  asymptotes.  The 
vertical  asymptotes  are  given  by  eq.  40  and  the  slant  asymptotes  by  eq.  42.  The  eigen¬ 
value  equation,  eq.  28,  is  solved  by  iteration,  using  eq.  29  and  iterating  on  x  for  the  A 
modes  and  using  eq.  30  and  iterating  on  y  for  the  B  modes.  The  iteration  is  started  at 
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the  horizontal  asymptotes  for  large  values  of  minus  x  for  the  B  modes  and  at  the 
vertical  asymptotes  for  large  values  of  minus  y  for  the  A  modes. 


Figure  14.  Eigenvalues  for  S  =  0.5, 
p  =  1.2,  free-rigid  boundary. 


Consider  now  the  A  modes  with  y  >  0.  These  are  modes  with  phase  velocities 
larger  than  the  barrier  sound  speed.  These  modes  correspond  to  rays  reflecting  off 
both  the  upper  and  lower  boundaries.  There  are  no  coupling  effects  for  these  modes. 
We  observe  that  here,  the  eigenvalues  appear  to  be  about  evenly  spaced  as  is  the  case 
with  single  ducts  with  no  coupling.  The  eigenvalues  for  the  B  modes  also  show  even 
spacing  characteristics  of  modes  in  a  single  refractive  duct  with  boundaries. 

Our  main  interest  here  is  in  the  A  modes  with  y  <  0.  These  modes  correspond  to 
modes  trapped  in  either  the  upper  or  lower  ducts.  Modes  trapped  in  the  upper  duct 
reflect  from  the  surface  and  form  a  refractive  nadir  in  the  upper  duct.  These  modes 
may  be  identified  in  figure  14  as  those  that  follow  the  vertical  asymptotes.  Modes 
trapped  in  the  lower  duct  reflect  from  the  rigid  bottom  and  form  a  refractive  apex  in 
the  lower  duct.  These  modes  may  be  identified  in  figure  14  as  those  that  follow  the 
slant  asymptotes. 

Observe  that  in  contrast  to  the  B  modes  or  the  A  modes  with  y  >  0,  the  A  modes 
for  y  <  0  have  a  highly  variable  spacing.  This  variability  represents  the  main  reason 
why  the  determination  of  modes  in  a  double  duct  is  so  complicated  as  compared  to 
that  for  regularly  spaced  modes  in  a  single-duct  configuration. 
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The  intersection  of  the  vertical  with  the  slant  asymptotes  correspond  to  the  critical 
frequencies.  At  or  near  these  points  the  adjacent  modes  are  transitioning  from  being 
trapped  in  the  upper  duct  to  being  trapped  in  the  lower  duct  and  vice  versa.  The  maxi¬ 
mum  coupling  between  upper  and  lower  ducts  occurs  near  these  points.  For  ease  of 
identification,  these  intersections  have  been  identified  by  Roman  numerals. 

We  note  that  the  eigenvalue  curves  in  figure  14  never  cross  each  other.  They  come 
very  close  to  each  other  near  the  intersection  points.  In  general,  the  more  negative  the 
y  of  the  intersection  point,  the  closer  the  adjoining  modes  come  together.  Not  all  inter¬ 
section  points  lead  to  critical  frequencies.  Critical  frequencies  do  not  occur  when  the  y 
of  the  intersection  point  is  positive.  Examples  in  figure  14  are  the  intersection  of  slant 
asymptote  two  with  vertical  asymptotes  one  or  two. 

A  higher  order  mode  may  have  many  critical  frequencies  associated  with  it.  Con¬ 
sider  for  example  mode  4  in  figure  14.  At  large  negative  y,  this  is  on  upper  duct 
mode.  At  critical  point  V,  it  transitions  to  a  lower  duct  mode.  At  IV,  it  transitions  back 
to  an  upper  duct  mode.  At  IE,  it  transitions  again  to  a  lower  duct  mode. 

We  note  that  for  sufficiently  large  negative  y  (high  frequency)  all  modes  start  out 
as  lower  duct  modes,  i.e.,  start  out  on  the  vertical  asymptotes.  Since  S  is  positive,  C, 
in  figure  1  is  less  than  C3.  At  sufficiently  high  frequency,  the  phase  velocity  for  all 
modes  goes  to  C,  since  it  is  less  than  C3.  Thus,  at  sufficiently  high  frequencies,  any 
mode  is  an  upper  duct  mode. 

The  eigenvalue  solutions  of  figure  14  are  not  exclusive  to  underwater  acoustics. 
They  represent  the  solution  to  all  equivalent  physical  problems  with  the  same  boundary 
conditions  and  canonical  parameters.  This  plot  could  apply  for  example  to  air  acous¬ 
tics,  seismic  propagation,  electromagnetics,  solid  state  physics,  and  perhaps  to  quantum 
mechanics.  The  particular  physical  application  only  comes  into  play  when  we  convert 
canonical  functions  into  the  physical  solutions. 

Figure  15  presents  the  sound-speed  profile  used  in  the  evaluation  of  the  various 
physical  quantities  of  this  section.  The  layer  1  parameters  are  the  same  as  those  dis¬ 
cussed  for  the  left  panel  of  figure  3.  These  parameters  and  the  values  of  p  and  S 
completely  determine  the  sound-speed  profile. 

Figure  16  presents  the  canonical  frequency  versus  x  corresponding  to  the  eigen¬ 
values  of  figure  14.  The  conversion  to  acoustical  frequency  is  for  the  profile  parame¬ 
ters  of  figure  15.  Later  on,  eigenfunctions  will  be  presented  for  the  critical  frequencies 
I,  n,  and  III  as  well  as  for  the  frequencies  of  500  and  900  Hz.  The  vertical  portions  of 
these  curves  correspond  to  modes  trapped  in  the  upper  duct.  Thus,  modes  1  and  3  at 
500  Hz  and  modes  1,  2,  and  4  at  900  Hz  are  trapped  in  the  upper  duct.  Portions  of 
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Figure  15.  Sound-speed  profile  used 
in  the  analysis  of  the  A  modes  of 
figure  14. 
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Figure  16.  Canonical  or  acoustic  fre¬ 
quencies  for  the  A  modes  of  figures 
14  or  15. 


the  curves  slanting  from  left  to  right  correspond  to  modes  trapped  in  the  lower  duct. 
Examples  are  modes  2  and  4  at  500  Hz  and  modes  3  and  5  at  900  Hz.  Mode  5  at  500 
Hz  as  well  as  other  modes  below  the  knee  in  the  curves  generally  correspond  to  rays 
reflecting  from  both  interfaces. 

Figure  17  present:;  the  phase  velocity  versus  frequency.  The  modes  with  phase 
velocity  less  than  the  dotted  line  at  1482  m/s  correspond  to  modes  trapped  in  the 
upper  or  lower  ducts.  The  phase  velocities  exhibit  the  characteristic  pinching  together 
at  the  critical  frequencies  I,  II,  and  ID.  At  sufficiently  high  frequencies,  all  modes  will 
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be  trapped  in  the  upper  duct  and  will  go  to  1480  m/s  as  a  limit.  Modes  trapped  in  the 
lower  duct  tend  to  a  limit  toward  1481  m/s,  which  is  the  bottom  sound  speed.  Exam¬ 
ples  are  mode  2  between  frequency  I  and  II  and  mode  3  beyond  frequency  II. 


Figure  17.  Phase  velocities  for  the  A 
modes  of  figure  15. 


Figure  18  presents  the  canonical  group  velocity  factor  of  eq.  58.  It  appears  rather 
messy,  but  is  just  what  is  required  to  generate  the  group  velocity  curves  about  to  be 
presented.  The  segments  with  near-constant  value  of  0.5  near  the  top  of  the  plot  are 
associated  with  the  slant  asymptotes  of  figure  14  having  a  constant  slope  of  -1. 

Figure  19  presents  group  velocity  versus  frequency.  At  the  critical  frequencies,  the 
group  velocities  of  adjacent  modes  cross  each  other.  Group  velocity  is  a  much  more 
sensitive  indicator  of  the  nature  of  the  modes  than  is  the  phase  velocity.  Modes 
trapped  in  the  upper  duct  tend  toward  1480  m/s,  the  surface  sound  speed.  Modes 
trapped  in  the  lower  duct  have  group  velocities  exceeding  1481  m/s,  the  bottom  sound 
speed.  The  rapid  drop  in  group  velocities  at  the  lower  frequency  end  of  the  plot  repre¬ 
sents  modes  corresponding  to  rays  reflecting  from  both  boundaries.  We  note  that  fig¬ 
ure  19  exhibits  many  points  of  Airy  phase  (extremes  in  the  group  velocity).  The  sig¬ 
nificance  of  the  Airy  phase  as  related  to  double-duct  propagation  is  discussed  on  pages 
41  and  42  of  reference  4. 
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Figure  18.  Group  velocity  factor 
for  the  A  modes  of  figure  14. 
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Figure  19.  Group  velocities  for 
the  A  modes  of  figure  15. 
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Figure  20  presents  the  mode  normalization  coefficient.  The  left  scale  presents  the 
canonical  normalization  factor.  The  right  scale  presents  the  normalization  coefficients 
as  obtained  by  scaling  with  the  constant  J  of  eq.  86.  The  bottom  scale  presents  the 
canonical  frequency.  The  top  scale  presents  the  acoustic  frequency  obtained  with  the 
use  of  the  scale  constant  K  of  eq.  55.  Observe  that  the  curve  for  mode  2  goes  off  the 
top  scale  at  a  frequency  of  about  300  Hz  and  returns  at  a  frequency  of  about  800  Hz. 
In  this  interval,  this  curve  forms  a  relative  maximum  with  Dn  =  8.9  x  105.  This  interval 
corresponds  to  the  interval  in  figure  14  where  the  second  mode  is  trapped  in  the  lower 
duct.  In  our  approach,  we  have  chosen  the  unnormalized  depth  functions  to  automati¬ 
cally  satisfy  the  boundary  condition  at  the  surface.  This  choice  results  in  huge  values 
for  the  normalization  coefficient  for  modes  when  they  are  trapped  in  the  lower  duct. 
This  does  not  present  a  problem  with  evaluating  mode  amplitudes,  but  only  with  the 
presentation  of  Dn  versus  frequency  on  a  single  scale. 

The  curves  go  off  scale  for  each  mode  when  it  is  trapped  in  the  lower  duct.  For 
example,  the  curve  for  mode  3  goes  off  scale  near  800  Hz.  Similarly,  the  curve  for 
mode  5  is  off  scale  between  800  and  1250  Hz.  Aside  from  these  relative  maxima  the 
curves  have  a  consistent  behavior.  The  normalization  curves  for  all  modes  tend  to  zero 
at  low  frequencies.  At  high  frequencies,  all  modes  trapped  in  the  upper  duct  appear  to 
merge  together  in  a  common  curve  with  a  small  negative  slope. 

The  curves  for  adjacent  modes  cross  each  other  at  the  critical  frequencies.  Thus, 
in  figure  20  modes  1  and  2  cross  at  frequency  I  and  modes  4  and  5  cross  at  fre¬ 
quency  in. 

Figure  21  presents  the  eigenfunctions  at  500  Hz  for  modes  1  to  5.  These  are  the 
normalized  mode  amplitude  functions  plotted  as  a  function  of  receiver  depth  from  the 
surface  to  the  bottom  at  186  m.  The  source  depth  is  45.4  m  for  modes  1,  3,  and  5. 
Corresponding  values  for  modes  2  and  4  are  185.9  and  140  m.  The  dotted  horizontal 
line  represents  the  interface  between  ducts  at  99.8  m.  The  dotted  vertical  line  is  the 
zero  axis  for  the  amplitudes. 

All  mode  amplitudes  shown  here  and  in  similar  figures  to  follow  in  this  section 
have  certain  features  in  common.  Because  of  the  boundary  conditions,  all  mode  ampli¬ 
tudes  have  a  node  at  the  surface  and  a  relative  maximum  at  the  bottom  boundary. 
Moreover,  the  modes  of  order  n  all  have  n-1  zero  crossings.  These  features  are  usually 
apparent  in  Figs.  21  and  others  to  follow.  In  a  few  cases,  the  amplitudes  are  too  small 
to  see  their  features.  In  these  cases,  they  have  been  verified  by  examining  computer 
printouts  of  the  numerical  values. 
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Figure  20.  Canonical  or  customary 
eigenvalue  normalization  coefficients 
for  the  A  modes  of  figures  14  or  15. 
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Figure  21.  Eigenfunctions  at 
500  Hz  for  the  A  modes  of  fig¬ 
ure  3. 
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Clearly,  Figure  21  shows  that  modes  1  and  3  are  trapped  in  the  upper  duct, 
whereas  modes  2  and  4  are  trapped  in  the  lower  duct.  (We  note  that  these  results  were 
predicted  in  the  discussion  of  Figs.  16  and  14.)  In  this  example,  mode  5  has  compara¬ 
ble  amplitudes  in  both  ducts.  The  y  eigenvalue  for  mode  5  is  positive.  Thus,  this  mode 
corresponds  to  rays  reflecting  from  both  interfaces.  This  is  of  no  particular  interest  in 
this  paper.  Our  chief  interest  is  in  comparable  amplitudes  in  the  two  ducts  resulting 
from  barrier  penetration  at  the  critical  frequencies  by  the  modes  normally  trapped  in 
the  individual  ducts. 

Figure  22  presents  the  eigenfunctions  at  900  Hz  for  modes  2  to  5.  We  have  not 
treated  mode  1  highly  trapped  in  the  upper  duct  because  of  a  problem  with  inadequate 
precision  in  the  numerical  evaluation.  The  source  depths  for  modes  2  to  5  were  respec¬ 
tively  45.4,  185.9,  60.0,  and  185.9  m.  These  depths  were  chosen  to  enhance  the  mode 
amplitudes  so  that  the  features  of  the  lower  amplitudes  could  be  shown.  Clearly, 
modes  2  and  4  are  upper  duct  modes,  whereas  modes  3  and  5  are  lower  duct  modes 
as  predicted  in  the  discussion  of  figures  16  and  14.  We  note  in  all  cases  of  trapped 
modes,  the  amplitudes  are  large  in  one  duct  but  not  in  the  other.  Thus,  the  response 
for  a  source  in  one  duct  and  receiver  in  the  other  duct  is  poor  for  trapped  modes. 
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Figure  22.  Eigenfunctions  at  900 
Hz  for  the  A  modes  of  figure  3. 


We  will  now  look  at  the  mode  amplitudes  at  the  critical  frequencies.  These  are 
frequencies  where  mode  pairs  are  transitioning  from  being  trapped  in  one  duct  to 
being  trapped  in  the  other  duct.  For  these  frequencies  the  amplitudes  are  comparable 
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for  the  two  ducts,  even  though  the  mode  phase  velocity  is  well  below  the  barrier  sound 
speed. 

Figure  23  presents  the  eigenfunctions  for  modes  1  and  2  at  269  Hz,  as  associated 
with  critical  frequency  I.  The  source  depth  is  39.4  m.  The  frequency  of  269  Hz  was 
determined  empirically  so  as  to  match  the  amplitudes  in  the  upper  duct.  The  dashed 
curve  is  the  result  of  folding  over  the  negative  portion  of  mode  2.  In  this  example, 
there  is  a  reasonably  good  match  between  the  amplitudes  of  modes  1  and  2. 
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Figure  23.  Eigenfunctions  at 
269  Hz  (critical  frequency  I). 


We  note  that  in  figures  21  and  22,  the  mode  number  is  given  by  the  number  of 
antinodes.  However,  in  figure  23,  mode  1  has  two  antinodes.  In  this  example  and 
others  to  follow,  in  the  regime  about  a  critical  frequency,  the  lower  order  mode  will 
have  an  extra  antinode.  With  a  knowledge  of  the  boundary  conditions,  the  mode  num¬ 
ber  may  often  be  determined  from  the  number  of  nodes  or  antinodes.  However,  near 
critical  frequencies,  an  extra  antinode  is  formed.  Thus,  the  only  fail-safe  method  is  for 
the  mode  number  to  be  always  one  more  than  the  number  of  zero  crossings. 

The  cross-duct  propagation  between  mode  1  and  2  is  enhanced  by  placing  the 
source  at  the  antinode  in  the  upper  layer  and  the  receiver  at  the  bottom  of  the  lower 
layer.  However,  in  this  configuration,  modes  1  and  2  are  of  opposite  phase  and  will 
cancel  each  other  at  zero  range.  The  modes  will  periodically  reinforce  each  other  at 
odd  half  multiples  of  the  mode  interference  range  and  will  cancel  each  other  at  even 
half  multiples.  This  is  true  of  the  next  two  examples  as  well.  The  difference  between 
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the  mode  phase  velocities  for  this  example  is  0.08569  m/s.  At  this  frequency  this 
translates  to  an  interference  range  of  95.2  km  as  given  by  eq.  117.  Figures  4-7  of 
reference  4  present  examples  of  the  interference  patterns  assoc:. ■  :  a  critical 
frequency. 

Figure  24  presents  the  eigenfunctions  for  modes  4  and  5  at  758  Hz,  as  associated 
with  critical  frequency  HI.  The  source  depth  is  35.4  m.  The  da  >iied  line  is  the  result  of 
folding  over  the  mode  5  curve  below  the  fourth  node.  Here,  the  match  between  the 
amplitudes  of  modes  4  and  5  is  somewhat  better  than  the  match  of  figure  23. 
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Figure  24.  Eigenfunctions  at 
758  Hz  (critical  frequency  ID). 


Optimum  cross-duct  propagation  is  obtained  by  placing  the  source  at  the  third 
antinode  in  the  upper  duct  and  the  receiver  at  the  fourth  antinode  of  mode  5  in  the 
lower  duct.  Again,  the  modes  are  out  of  phase  below  the  duct.  The  modes  will  again 
reinforce  or  cancel  each  other.  The  difference  between  the  mode  phase  velocities  is 
0.02026  m/s,  corresponding  to  an  interference  range  of  143.0  km. 

Figure  25  presents  the  eigenvalues  for  modes  2  and  3  near  799  Hz,  as  associated 
with  critical  frequency  II.  The  source  depth  is  45.4  m.  The  dashed  line  is  the  result  of 
folding  over  the  mode  3  curve  below  the  third  node.  Here,  there  is  almost  a  perfect 
match  between  the  eigenvalues  because  the  phase  velocities  lie  so  close  to  each  other. 
Orthogonality  between  the  modes  is  preserved  by  the  fact  that  the  modes  have  opposite 
phase  below  the  third  node  of  mode  3. 
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Figure  25.  Eigenfunctions  at 
798.921  Hz  (critical  frequency  13). 


Optimum  cross-duct  propagation  is  obtained  by  placing  the  source  at  the  second 
antinode  in  the  upper  duct  and  the  receiver  at  the  bottom  of  the  lower  layer.  The 
difference  between  mode  phase  velocities  is  only  0.00012  m/s.  This  results  in  a  whop¬ 
ping  interference  range  of  2.3  x  104  km.  Such  closely  coupled  modes  have  been 
assessed  in  reference  4.  When  properly  summed  up,  the  two  modes  do  not  make  a 
significant  contribution  to  the  propagation  between  the  upper  and  lower  ducts. 

C.  FREE-FREE  BOUNDARY,  p  =  1,  S  =  0 

This  section  as  well  as  the  next  three  sections  treat  the  case  of  p  =  0,  which  occurs 
when  the  sound  speed  is  the  sqme  at  the  upper  and  lower  boundaries.  This  case  is  of 
particular  interest  because  of  past  problems  with  the  phase  integral  approach  discussed 
in  connection  with  eq.  68  and  its  free-free  and  rigid-rigid  boundary  counterparts.  For 
the  free-free  boundary  eq.  68  becomes 

P  =  3m/3n  • 

This  section  treats  the  case  of  p  =  1,  i.e.,  m  =  n.  This  is  the  case  of  all  two-layer 
profiles  symmetric  about  the  interface  between  layers  1  and  2  and  having  both  a  free- 
surface  and  free-bottom  condition.  Figure  26  presents  the  canonical  eigenvalues  for  this 
case.  The  method  of  designating  modes  and  asymptotes  has  been  changed  from  that  of 
figure  2  and  most  other  figures  and  the  number  of  modes  shown  has  been  increased 
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from  1  to  5  to  1  to  10.  These  changes  are  made  to  facilitate  the  analysis  of  section  3E 
for  variable  p.  Rather  than  have  a  distinctive  curve  for  modes  1  to  5  as  in  most  fig¬ 
ures,  figure  26  displays  the  curves  for  odd  modes  by  a  long  dash  separated  by  a 
shorter  dash.  The  even  modes  are  shown  by  the  curves  with  short  dashes. 


|  T - 1 - 1  I  I  1  I - 1 
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Figure  26.  Eigenvalues  for  S  =  0,  p  =1,  free-free 
boundary. 


There  are  three  major  differences  from  the  free-rigid  counterpart  of  figure  2.  The 
most  significant  difference  is  that  the  two  sets  of  vertical  asymptotes  of  figure  3  col¬ 
lapse  to  one  set  in  figure  26,  given  by  the  solid  vertical  lines.  Thus,  the  modes  do  not 
go  to  separate  asymptotes  as  in  figure  2,  but  pairs  of  modes  go  to  a  common  asymp¬ 
tote.  The  odd  modes  approach  the  asymptotes  from  the  left,  whereas  the  even  modes 
approach  from  the  right.  Thus,  the  eigenvalues  of  modes  1  and  2,  modes  3  and  4,  etc 
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have  common  values  of  x  =  ai,  a2,  etc.  and  y  =  -  oo.  Physically,  this  means  the  eigen¬ 
values  converge  at  infinite  frequency. 

Another  significant  difference  is  that  in  figure  26,  the  line  y  =  x  is  a  solution  to  the 
eigenvalue  equation,  whereas  in  figure  2,  this  line  is  not.  Note  that  for  values  of  x  and 
y  greater  than  about  1,  the  eigenvalue  curves  have  almost  uniform  spacing.  For  exam¬ 
ple,  the  spacing  between  modes  1A  and  2A  is  comparable  to  the  spacing  between 
mode  1A  and  y  =  x  or  between  modes  IB  and  2B  or  between  mode  IB  and  y  =x.  In 
contrast,  in  figure  2  the  line  y  =x  is  not  a  solution  and  the  spacing  between  modes  1A 
and  IB  is  comparable  to  that  between  modes  1A  and  2A  and  IB  and  2B.  Thus,  the 
status  of  the  line  y  =  x  explains  the  general  location  of  the  eigenvalues  for  mode  1. 
This  explanation  is  not  available  in  the  conventional  eigenvalue  approach,  where  the  A 
and  B  modes  are  treated  as  independent  problems  and  the  solution  y  =  x  has  no  coun¬ 
terpart  at  all. 

A  third  difference  is  the  fact  that  the  curves  of  figure  26  go  through  two  sets  of 
lattice  points.  One  of  these  sets  is  given  by  the  intersection  of  the  solid  grid  lines,  i.e., 
the  solutions  given  by  eq.  47.  The  other  set  is  given  by  the  intersection  of  the  dotted 
grid  lines,  i.e.,  the  solutions  given  by  eq.  48.  The  points  at  x  =  -  <»  or  y  =  -  corre¬ 
sponding  to  the  asymptotes,  are  lattice  points  for  the  solid  grid,  but  not  for  the  dotted 
grid.  Notice  that  each  lattice  point  lies  on  one  of  the  modal  curves  or  on  the  line 
y  =  x.  Only  one  curve  goes  through  each  lattice  point  except  for  the  solutions  at 

y  =  -  oo. 

There  are  no  such  lattice  points  in  figure  2.  This  difference  has  no  particular  sig¬ 
nificance.  In  the  case  of  the  one-layer  ducts  of  reference  2,  all  boundary  conditions 
lead  to  eigenvalue  curves  with  this  lattice  point  property.  In  the  case  of  two-layer  ducts, 
only  the  case  of  figure  26  and  the  corresponding  case  with  rigid-rigid  boundary  condi¬ 
tions  (presented  in  section  3D)  have  lattice  points. 

We  now  proceed  to  evaluate  the  case  A  modes  for  the  standard  profile  parameters 
used  in  Section  2A.  The  sound-speed  profile  is  identical  to  case  A  of  figure  3  except 
with  a  free  rather  than  rigid  bottom.  Figures  27  to  32  are  the  free-free  counterparts  of 
figures  4  to  9.  These  modes  are  designated  by  the  general  identification  scheme  and 
do  not  follow  the  special  scheme  of  figure  26.  The  curves  of  figures  27,  28,  30,  and 
31  merge  together  in  pairs  at  high  frequencies  in  contrast  to  their  free-rigid  counter¬ 
parts  in  which  all  modes  are  well  separated.  The  curves  of  figure  29  connect  together 
in  pairs  in  contrast  to  figure  6  where  they  remain  distinct.  We  note  that  Dn  in  fig¬ 
ure  31  does  not  become  very  large  as  is  the  case  for  the  odd  modes  of  figure  8.  This 
problem  disappears  when  the  boundary  conditions  are  the  same  at  the  surface  as  at 
the  bottom. 
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Figure  31.  Canonical  or  customary  eigenvalue  nor¬ 
malization  coefficients  for  the  A  modes  of  figure  26. 

Figure  32  presents  the  eigenfunctions  for  modes  1  to  5  at  a  frequency  of  500  Hz. 
The  source  depth  is  25,  25,  60,  135,  and  100  for  modes  1  to  5  respectively.  Again,  the 
mode  number  is  one  more  than  the  number  of  zero  crossings.  Here,  the  number  of 
modes  is  n  +  1  where  n  is  the  mode  number.  The  number  of  antinodes  is  equal  to  the 
mode  number  for  even  modes  and  is  equal  to  n  +  1  for  odd  modes.  We  note  that  the 
eigenfunctions  are  symmetric  about  the  interface  between  layers  for  odd  modes  and 
antisymmetric  for  even  modes. 

As  can  be  seen  in  figure  32,  there  is  a  strong  resemblance  between  the  eigenfunc¬ 
tions  for  mode  1  and  2.  As  discussed  in  a  later  section,  the  coupling  between  the 
upper  and  lower  duct  at  this  frequency  appears  to  be  in  line  with  other  examples  of 
critical  frequency.  While  this  example  has  coupling  features  in  common  with  those  of 
critical  frequencies  for  S  *  0,  there  are  significant  differences.  When  S  *  0,  there  are 
critical  frequencies  for  a  given  mode  when  paired  with  all  other  modes.  In  this  case  of 
S  =  0,  each  odd  mode  is  paired  with  the  next  higher  even  mode.  There  are  no  discrete 
critical  frequencies  as  for  S  *  0.  For  each  pair  of  odd-even  modes  the  modes  become 
coupled  at  sufficiently  high  frequency  with  the  coupling  increasing  as  the  frequency  is 
increased.  At  infinite  frequency,  the  eigenvalues  become  identical  for  each  pair  of 
odd-even  modes.  An  important  feature  is  that  these  pairs  have  different  eigenvalues  for 
all  finite  frequencies.  This  ic  in  marked  contrast  to  the  phase-integral  method  which 
predicts  pairs  with  identical  eigenvalues  (degenerate  modes)  for  all  frequencies. 
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Figure  32.  Eigenfunctions  at  500  Hz  for  the  A  modes  of  figure  26. 


D.  RIGID-RIGID  BOUNDARY,  p  =  1,  S  =  0 

This  case  is  the  rigid-rigid  counterpart  of  the  profile  of  section  3C  and  has  many 
similar  features.  Figure  33  presents  the  eigenvalues.  The  special  format  is  that  dis¬ 
cussed  for  figure  26.  In  this  case,  the  common  pair  of  vertical  asymptotes  are  given  by 
x  =  a',,  a'2,  etc.,  rather  than  the  a,,  a2,  etc.  of  figure  26.  The  line  y  =  x  is  a  solution 
as  in  figure  26.  The  curves  of  figure  33  go  through  the  lattices  points  of  eq.  49  (solid 
lines)  or  50  (dotted  lines)  rather  than  those  of  eqs.  47  or  48  in  figure  26. 
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Figure  33.  Eigenvalues  for  S  =  0,  p  =  1,  rigid-rigid  boundary. 
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The  most  significant  difference  between  figure  33  and  figure  26  is  the  peculiar 
behavior  of  mode  1  in  figure  33.  Here,  mode  1A  connects  with  mode  IB  at  x  =  0, 
y  =  0.  This  behavior  occurs  for  all  profiles  with  rigid-rigid  boundary  conditions.  This 
case  is  analyzed  in  detail  in  the  discussion  of  figure  4  of  reference  3,  which  treats  the 
single  layer  case.  The  chief  characteristic  is  that  mode  1A  or  IB  for  the  rigid-rigid 
boundary  condition  only  exists  for  refracted  rays.  This  condition  does  not  exist  for  rays 
reflecting  from  both  boundaries. 

Figure  34  presents  the  canonical  or  acoustic  frequency  plot.  The  distinctive  differ¬ 
ence  from  figure  27  is  that  mode  1  in  figure  34  has  no  lower  cutoff  frequency  as  do 
the  other  modes.  The  chief  difference  between  the  phase  velocities  of  figure  35  and 
those  of  figure  28  is  that  the  phase  velocity  of  mode  1  in  figure  35  does  not  go  to 
infinity  for  low  frequencies  by  approaches  a  constant  value  as  f  goes  to  zero.  This 
constant  value  is  given  by  eq.  36  of  reference  3  and  is  the  square  root  of  the  harmonic 
mean  of  the  square  of  C,  and  Q.  This  is  very  slightly  less  than  the  mean  of  C,  and 
C2,  which  in  this  case  is  1481  m/s. 


Figure  34.  Canonical  or  acoustic  frequencies 
for  the  A  modes  of  figure  33. 
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Figure  35.  Phase  velocities  for  the  A 
modes  of  figure  33. 


Figure  36  for  the  canonical  group  velocity  factor  is  similar  to  figure  29  with  the 
exception  of  mode  1.  It  does  not  go  to  large  negative  values  for  small  x  but  goes  to 
0.5  as  x  goes  to  zero.  The  analyses  of  reference  3  notes  that  the  slope  of  mode  1  at 
the  origin  in  figure  33  is  -1.  This  makes  the  group  velocity  factor  equal  to  2*1. 

Figure  37  for  the  group  velocity  is  similar  to  figure  30  with  the  exception  of  mode 
1.  As  f  goes  to  zero  the  group  velocity  does  not  go  to  zero  as  for  other  modes,  but 
goes  to  a  constant  that  is  the  same  as  the  limit  of  phase  velocity  of  figure  35. 

Figure  38  for  the  eigenvalue  normalization  coefficient  is  quite  different  from  fig¬ 
ure  31.  In  figure  31  the  coefficient  goes  to  zero  at  low  frequencies  and  that  for  all 
modes  merges  at  the  higher  frequencies.  In  figure  38,  mode  1  goes  to  one  constant  at 
zero  frequency  while  all  other  modes  go  to  another  constant.  This  same  behavior  is 
observed  for  the  A  modes  in  figure  22  of  reference  3.  Indeed  reference  3  demonstrates 
that  the  limits  for  mode  1  are  "ir2  =  0.101  and  2tt‘2  =  0.202.  In  figure  38  the  modes 
merge  in  pairs  at  higher  frequencies.  However,  the  pairs  do  not  also  merge  as  in  fig¬ 
ure  31.  This  is  similar  to  the  case  of  figure  22  of  reference  3  in  which  none  of  the  A 
modes  have  merged,  at  least  for  the  frequency  range  shown. 
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Figure  38.  Canonical  or  customary  eigenvalue 
normalization  coefficient  for  the  A  modes  of 
figure  33. 


Figure  39  presents  the  eigenfunctions  for  modes  1  to  5  at  a  frequency  of  400  Hz. 
The  source  depth  is  5,  5,  50,  50,  and  100  for  modes  1  to  5  respectively.  Again,  the 
mode  number  is  one  more  than  the  number  of  zero  crossing.  The  number  of  nodes  is 
one  less  than  the  mode  number.  The  number  of  antinodes  is  n  +  1  for  odd  modes  and 
is  equal  to  n  for  even  modes.  As  was  the  case  of  figure  32,  the  eigenfunction  are 
symmetric  about  the  interface  between  layers  for  odd  modes  and  antisymmetric  for 
even  modes. 

The  general  characteristics  are  essentially  the  same  as  those  discussed  in  the  last 
paragraph  of  section  3C. 


E.  EIGENVALUES  FOR  FREE-FREE  BOUNDARY,  VARIABLE  p,  S  =  0 

In  this  section,  we  determine  the  behavior  of  the  eigenvalue  curves  for  S  =  0  as  we 
let  p  decrease  or  increase  from  p  =  1.  We  chose  the  free-free  boundary  condition  for 
illustrative  purposes  as  once  this  is  available,  the  situation  for  other  boundary  condi¬ 
tions  can  be  easily  inferred. 
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Figure  39.  Eigenfunctions  at  400  Hz  for  the  A  modes  of  figure  33. 
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Figure  40  presents  the  eigenvalues  for  p  =  0.99.  The  format  is  similar  to  that  of 
figure  26.  The  dotted  grid  lines  of  figure  26  are  deleted  because  there  are  no  lattice 
point  solutions  for  p  *  1.  There  are  two  sets  of  vertical  asymptotes  in  figure  40  as 
compared  to  the  one  set  of  figure  26.  The  solid  vertical  asymptotes  are  the  same  as 
those  of  figure  25,  i.e.,  x  =  an.  The  eigenvalues  for  the  odd-order  modes  of  figure  40 
are  almost  the  same  as  those  of  figure  26  because  the  vertical  asymptotes  for  the  odd 
modes  remains  unchanged  for  different  p.  The  vertical  asymptotes  for  the  even-order 
modes  of  figure  40  are  shown  by  the  dotted  vertical  lines  as  given  by  x  =  an/p2.  The 
displacement  between  the  sets  of  asymptotes  for  odd  and  even  modes  is  given  by 
Ax  =  an(l  -  p'2).  Thus,  for  a  given  value  of  p2  this  displacement  increases  with 
increasing  mode  number  since  an  increases  with  mode  number.  If  p  is  decreased  fur¬ 
ther  below  1,  the  even  modes  will  march  further  to  the  right  than  illustrated  in  fig¬ 
ure  40. 
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Figure  40.  Eigenvalues  for  S  =  0,  p  =  0.99,  free-free 
boundary. 
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Next,  we  examined  a  value  of  p  =  1.01.  On  the  basis  of  figure  40  and  other 
smaller  values  of  p  we  assumed  that  the  odd  modes  would  still  remain  fixed  and  that 
the  even  modes  would  march  to  the  left.  We  were  looking  forward  to  the  possibility  of 
the  even  mode  eigenvalue  crossing  the  odd  mode  eigenvalues  and  producing  degenerate 
eigenvalues.  However,  as  shown  in  figure  41,  this  does  not  occur.  In  figure  41,  the 
roles  of  even  and  odd  modes  reversed  at  p  =  1  and  as  p  is  increased  above  one  the 
even  mode  asymptotes  remain  fixed  at  x  =  an,  whereas  the  odd  mode  asymptotes 
march  to  the  left  according  to  x  =  a^p2. 


O 


Figure  41.  Eigenvalues  for  S  =  0,  p  =  1.01,  free-free 
boundary. 
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Next,  consider  what  happens  for  even  larger  values  of  p.  Figure  42  presents  the 
eigenvalues  for  p  =  1.515.  Here,  the  dashed  lines  correspond  to  the  variable  symptotes, 
where  the  long  dashed  lines  correspond  to  the  fixed  asymptotes.  The  salient  feature 
here  is  that  the  fixed  vertical  asymptotes  of  x  =  a,  and  x  =  a2  are  occupied  by 
modes  3  and  8  respectively,  whereas  modes  1,  2,  4  to  7,  9,  and  10  have  asymptotes  at 
x  =  a^p2  for  n  =  1  to  8  respectively.  To  explain  the  results  of  figure  42,  we  first 
examine  where  the  asymptotes  of  mode  2n  +  1  in  figure  41  first  encounter  mode  2n  as 
p  is  increased  above  1.01.  This  solution  is  given  by  p2  =  an+i/an.  Thus,  modes  3  and  2 
will  have  a  common  asymptote  at  x  =  2.338  when  p  =  1.3222712. 
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Figure  42.  Eigenvalues  for  S  =  0,  p  =  1.515,  free-free 
boundary. 
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As  p  is  increased  above  1.01,  i.e.,  the  case  of  figure  41,  mode  2  remains  fixed  at 
2.338  and  mode  3  moves  to  the  left  until  p  =  1.3222712.  As  p  is  further  increased, 
mode  3  now  remains  fixed  at  2.338  while  mode  2  continues  to  move  to  the  left  as  in 
the  case  of  p  =  1.515  in  figure  42.  When  p2  =  5.521/2.338  or  p  =  1.537,  modes  3  and 
4  will  have  a  common  asymptote  at  2.338  and  as  p  increases  beyond  1.537,  mode  4 
will  remain  fixed  at  2.338  while  mode  3  moves  to  the  left. 

In  figure  42,  we  see  that  mode  8  corresponds  to  the  fixed  asymptote  at  4.088. 

Mode  8  took  over  this  asymptote  when  p2  =  a6/a2  or  p  =  1.486.  Mode  9  will  take  over 
this  asymptote  when  p2  =  a7/a2  or  p  =  1.567. 

A  similar  situation  occurs  when  p  is  decreased  further  from  the  value  of  0.990  in 
figure  40.  The  asymptotes  of  mode  2n  first  encounters  mode  2n  +  1  when  p2  =  an/an+]. 
Thus,  mode  2  will  have  a  common  asymptote  at  x  =  4.088  when  p2  =  2.338/4.088  or 
p  =  0.756.  This  is  the  reciprocal  of  the  value  of  p  for  which  mode  3  moves  to  mode  2 
in  figure  41.  The  salient  feature  is  that  the  mode  eigenvalues  never  cross  each  other. 
For  special  values  of  p  they  have  a  common  asymptote,  i.e.,  they  touch  each  other  at 
y  =  -  oo  corresponding  to  infinite  frequency.  As  the  value  of  p  is  increased  from  these 
special  values,  one  asymptote  remains  fixed  while  the  other  moves  to  the  left.  As  the 
value  of  p  is  decreased  from  these  special  values,  one  asymptote  remains  fixed  while 
the  other  moves  to  the  right. 

A  similar  type  of  behavior  occurs  for  the  rigid-rigid  boundary  condition.  The  only 
difference  is  that  the  asymptotes  are  associated  with  an'  rather  than  an. 

The  behavior  for  variable  p  for  the  free-rigid  boundary  conditions  is  similar  to  that 
of  the  free-free  and  rigid-rigid  boundary  conditions.  As  p  is  increased  from  the  special 
values  given  by  eq.  71,  one  asymptote  remains  fixed  while  the  other  moves  to  the 
right.  The  significant  difference  for  the  free-rigid  condition  is  that  p  =  1  is  not  a  spe¬ 
cial  value  as  for  the  free-free  and  rigid-rigid  conditions.  Indeed,  the  maximum  separa¬ 
tion  for  mode  eigenvalues  occur  for  p  =  1  as  illustrated  by  figure  2. 

It  is  of  interest  to  examine  what  happens  as  p  goes  to  its  limiting  values.  The  same 
behavior  applies  to  all  three  boundary  conditions.  We  have  no  interest  in  negative  p 
because  this  does  not  represent  a  double-duct  case.  As  p  decreases  from  1  to  zero  as  a 
limit,  the  variable  asymptotes  correspond  to  the  higher  order  modes  and  all  march  out 
toward  x  =  co,  leaving  the  fixed  asymptotes  to  the  lower  order  modes.  As  p  is 
increased  from  1  to  oo  as  a  limit,  the  variable  asymptotes  correspond  to  the  low-order 
modes  and  all  march  toward  x  =  0,  leaving  the  fixed  asymptotes  to  the  high-order 
modes. 

F.  FREE-RIGID  BOUNDARY,  p  =  0.6601015,  S  =  0. 

For  this  example,  p  satisfies  eq.  68  with  m  and  n  =  1.  This  is  a  case  where  modes 
1  and  2  have  a  common  vertical  asymptote  of  x  =  2.338. 
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This  example  was  prompted  by  an  earlier  investigation  of  a  positive-gradient  sur¬ 
face  duct  overlaying  a  symmetric  refractive  duct  where  the  surface  and  axial  sound 
speeds  were  the  same.  This  investigation  has  already  been  introduced  in  the  discussion 
of  eq.  69,  which  gives  the  relationship  between  gradients  designed  to  superimpose  the 
eigenvalues  for  mode  1  of  the  two  ducts.  The  unpublished  results  of  this  investigation 
displayed  an  unexpected  result.  At  400  Hz,  the  phase  velocities  for  modes  1  and  2 
differed  by  only  2.9  x  10*6  m/s.  However,  there  were  large  differences  in  the  eigen¬ 
functions.  The  mode  2  peak  in  the  upper  duct  had  45  dB  less  loss  than  the  mode  1 
peak.  Similarly,  the  mode  1  peak  in  the  lower  duct  had  47  dB  less  loss  than  the  mode 
2  peak. 

The  motivation  here  is  to  examine  a  case  where  modes  1  and  2  have  a  common 
asymptote.  The  cases  of  sections  3C  and  3D  have  this  property.  However,  modes  1 
and  2  strongly  resemble  each  other  in  figures  32  and  39.  Thus,  we  turn  to  the  case 
where  modes  1  and  2  have  a  common  asymptote  for  the  free-rigid  boundary. 

Figure  43  presents  the  canonical  eigenvalues.  Figure  44  gives  the  sound-speed  pro¬ 
file.  Figures  45  to  49  give  the  results  for  frequency,  phase  and  group  velocity,  and 
normalization  coefficients.  The  mode  1  and  2  results  of  figures  43,  45,  46,  47,  and  48 
are  in  general  qualitative  agreement  with  the  results  of  figures  26  to  30  respectively  for 
the  free-free  profile.  Moreover,  the  results  of  figures  46  and  48  are  in  general  agree¬ 
ment  with  the  phase  and  group  velocity  plots  of  reference  6,  i.e.,  the  results  of  the 
earlier  investigation  under  question. 

The  first  clue  that  there  is  a  significant  difference  is  figure  49.  Whereas  modes  1 
and  2  merge  at  high  frequencies  in  figure  31,  they  do  not  do  so  in  figure  49.  Indeed, 
computer  printouts  for  the  higher  frequencies  in  figure  49  indicate  that  the  mode  1 
values  continue  to  get  larger  with  increasing  frequency  until  computer  roundoff  prob¬ 
lems  set  in.  Thus,  we  are  not  certain  that  Dca  for  mode  1  remains  bounded  let  alone 
that  is  merges  with  the  relatively  small  values  for  mode  2. 

Figure  50  presents  the  eigenfunctions  for  modes  1  and  2  at  800  Hz.  The  source 
depth  is  25  m  .  The  difference  in  the  phase  velocity  of  modes  1  and  2  is  0.0061  m/s. 
We  see  that  at  this  frequency  mode  1  is  trapped  in  the  upper  duct  and  mode  2  is 
trapped  in  the  lower  duct. 

We  next  upped  the  frequency  to  bring  the  phase  velocities  differences  closer  with 
the  idea  of  bringing  the  eigenfunctions  into  better  agreement.  Figure  51  presents  the 
eigenfunction  for  1000  Hz.  Here  the  difference  in  the  phase  velocity  is  0.00012  m/s. 
However  the  amplitude  of  mode  1  in  the  upper  layer  is  even  smaller  than  in  figure  50. 
Detailed  computer  printouts  indicate  that  the  mode  1  antinode  in  the  upper  layer  has  a 
value  of  0.299  x  10*4.  This  is  62  dB  below  the  mode  2  antinode  in  the  upper  layer. 
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Figure  43.  Eigenvalues  for  S  =  0,  p  =  0.6601015, 
free-rigid  boundary. 
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Figure  44.  Sound-speed  profile  used  for 
the  analysis  of  the  A  modes  of  figure  43. 
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Figure  48.  Group  velocities  for  the  A  modes  of 
figure  43. 
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Figure  49.  Canonical  or  customary  eigenvalue 
normalization  coefficient  for  the  A  modes  of 
figure  43. 


Figure  50.  Eigenfunctions  at  800  Hz  for  modes 
1  and  2  for  the  A  modes  of  figure  43. 


Figure  51.  Counterpart  of  figure  50  at 
1000  Hz. 


Thus,  this  case  exhibits  the  same  type  of  behavior  encountered  in  the  earlier  investiga¬ 
tion,  i.e.,  very  similar  eigenvalues  but  radically  different  eigenfunctions  in  the  upper 
duct.  We  conclude  then  that  the  peculiar  results  of  the  earlier  investigation  are  correct 
and  not  due  to  computer  roundoff  or  other  errors  or  artifacts. 

Both  figures  50  and  51  exhibit  a  very  intriguing  result  in  that  the  amplitudes  for 
mode  1  and  2  are  almost  identical  in  the  lower  duct,  but  with  a  180-degree  change  in 
phase.  This  same  property  is  exhibited  in  all  our  examples  of  critical  frequencies  and 
in  figures  32  and  39  as  well.  Thus,  there  appears  to  be  some  “method  in  madness.” 
This  gives  us  confidence  that  the  results  of  figures  50  and  51  are  correct  and  not  the 
results  of  errors  or  artifacts. 

We  note  that  comparable  eigenfunctions  in  the  lower  duct  was  not  the  case  for  the 
earlier  investigation  that  obtained  radically  different  eigenfunctions  in  both  the  upper 
and  lower  ducts.  However,  this  result  should  not  be  surprising  because  of  the  different 
natures  of  the  lower  duct.  In  the  present  investigation,  the  lower  duct  is  a  negative 
gradient  layer  over  a  rigid  boundary.  In  the  earlier  investigation,  the  lower  duct  was  a 
refractive  duct  lying  over  a  negative  gradient  half  space. 

Since  we  have  no  physical  concept  that  explains  the  results  of  figures  50  and  51, 
we  turn  to  a  mathematical  analysis.  The  analysis  of  the  eigenvalues  itself  is  different 
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because  it  involves  not  only  an  analysis  of  Dn  but  also  of  the  unnormalized  eigen¬ 
values.  We  decided  to  determine  the  mathematical  circumstances,  which  lead  to  the 
divergence  of  Dca  between  modes  1  and  2  in  figure  49  and  compare  this  to  the  conver¬ 
gence  between  modes  1  and  2  in  figure  38.  The  rationale  for  this  comparison  is  the 
convergence  of  Dca  between  modes  1  and  2  in  figures  31  and  38  and  the  fact  that  Dca 
crosses  for  the  adjacent  modes  associated  with  the  critical  frequency  in  figure  19. 

These  results  are  also  associated  with  the  critical  frequency  in  figure  19.  These  results 
are  associated  with  eigenfunctions  which  strongly  resemble  each  other.  On  the  other 
hand,  the  divergence  in  figure  49  between  Dca  for  modes  1  and  2  is  clearly  related  to 
the  fact  that  the  eigenfunctions  of  modes  1  and  2  do  not  resemble  each  other. 

The  dominant  term  of  eq.  88  may  be  written  as 

Ai2(-  x)  [y  Bi2(-  y)  +  Bi'2(-y)]  . 

For  large  negative  y,  the  quantity  in  brackets  becomes  very  large.  However,  Ai(-x)  is 
near  zero  and  this  tends  to  keep  the  bracketed  term  in  check. 

In  the  case  of  1,000  Hz,  -x,  =  2.3377657  and  -x2  =  2.381100,  while  the  zero  of 
Ai(-x)  is  —  x0  =  2.3381074.  Thus,  -x2  -  (-x,)  =  0.0003443.  This  value  determines  the 
difference  in  phase  velocity,  i.e.,  the  closeness  of  the  eigenvalues.  However, 

-x0  -  (-x,)  =  0.0003417  and  -x0  -  (~x2)  =  -0.0000026.  Thus,  mode  2  lies  very  close  to 
the  zero,  whereas  mode  1  is  much  further  removed.  Indeed  [AK-x^/AiC-Xj,)]2  =  2.16  x 
106.  In  the  case  of  800  Hz,  this  ratio  =  5.98  x  103.  Thus,  although  x,  is  closer  to  x2 
and  to  x0  at  the  higher  frequency,  the  ifference  between  mode  1  and  2  in  the  normali¬ 
zation  term  continues  to  increase. 

Next,  we  do  a  similar  analysis  for  the  case  of  figure  38  at  400  Hz.  Here,  the  coef¬ 
ficient  of  the  bracketed  term  just  discussed  is  Ai'(-x).  Here,  -x,  =  1.018556,  -x2  = 
1.019010,  and  -x0  =  1.018793.  Thus  -x2  -  (-x,)  =  0.000454,  -x0  -  (x,)  =  0.000237, 
and  — x0  -  (-x2)  =  0.000217.  In  contrast  to  our  first  example,  modes  1  and  2  are  about 
equidistant  from  the  zero.  Here,  [Ai^-xO/Ai^-x^]2  =  1.19.  Thus,  not  surprisingly,  the 
modes  1  and  2  values  in  figure  38  are  close  together  at  400  Hz. 

Our  final  investigation  was  to  make  certain  that  the  type  of  behavior  exhibited  in 
figure  49  was  not  exclusive  to  the  free-rigid  boundary  condition.  Thus,  we  examined 
the  free-free  boundary  condition  for  S  =  0  and  p  =  1.322.  This  value  of  p  superimposes 
the  vertical  asymptotes  of  modes  2  and  3.  The  counterpart  of  figure  49  for  p  =  1.322 
demonstrated  that  the  Dca  for  mode  2  was  well  behaved  while  that  of  mode  3  diverged 
in  a  manner  similar  to  that  of  mode  1  in  figure  49.  We  conclude  then  that  the  behav¬ 
ior  encountered  here  for  the  free-rigid  boundary  condition  and  for  modes  1  and  2 
occurs  whenever  a  p  is  chosen  to  superimpose  vertical  asymptotes  for  any  two  modes 
and  for  any  boundary  condition. 
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G.  NEW  MEASURES  OF  MODE  COUPLING 

Reference  4  introduced  the  cross-duct  coupling  function  as  a  measure  of  mode 
coupling.  This  is  a  useful  measure  of  more  coupling  as  a  function  of  frequency  for  a 
given  profile  configuration.  However,  it  does  not  appear  to  be  as  useful  in  comparing 
different  profile  configurations  as  two  other  measures  of  mode  coupling  to  be  dis¬ 
cussed  here.  The  chief  purpose  for  introducing  these  other  measures  is  to  determine  if 
the  coupling  exhibited  in  figurs  32,  39,  50,  and  51  for  S  =  0  is  similar  to  or  different 
from  that  for  critical  frequencies  such  as  exhibited  in  figures  23,  24,  and  25. 

The  first  measure  of  coupling  that  comes  to  mind  is  the  difference  in  phase  velocity 
between  the  adjacent  modes,  because  in  general,  the  smaller  this  difference  the  closer 
the  modes  are  coupled.  However,  this  does  not  take  into  account  the  frequency 
dependence.  A  measure  which  takes  this  into  account  is  the  difference  in  wavelength 
or  An  =  ft>Cpn  Cpn+^ACm)'1  ■  However,  a  related  measure  is  the  mode  interference  range 
as  given  by  Rn+]  =  2irAA„.  Here,  the  closer  the  mode  coupling,  the  larger  the  interfer¬ 
ence  range. 

We  now  define  a  new  measure  of  mode  coupling.  For  our  examples  of  closely 
coupled  modes,  the  mode  amplituoe  of  the  lower  order  mode  formed  a  relative  mini¬ 
mum  whereas  the  higher  order  mode  had  a  zero  crossing.  We  observed  that  the  ampli¬ 
tude  at  the  relative  minimum  was  smaller  for  the  more  highly  coupled  modes.  Our  new 
measure  then  was  the  amplitude  of  this  relative  minimum  divided  by  the  largest  ampli¬ 
tude  for  this  mode.  We  note  that  this  measure  is  independent  of  source  depth,  because 
the  source  depth  factor  scales  both  amplitudes  by  the  same  amount. 

We  now  examine  the  correlation  between  our  new  measure  of  mode  coupling  and 
the  mode  interference  range.  This  correlation  is  shown  in  figure  52  on  a  log-log  scale. 
The  points  are  labeled  by  the  figure  number  corresponding  to  the  evaluation.  We  see 
that  the  points  for  figures  32  and  39  exhibit  the  same  general  relationship  as  that  of 
figures  23,  24,  and  25.  We  conclude  then  that  the  coupling  exhibited  in  figures  32 
and  39  is  essentially  the  same  as  that  for  critical  frequencies  as  exhibited  in  figures  23 
to  25.  The  evaluation  for  figure  50  is  given  by  the  cross.  The  corresponding  evaluation 
for  figure  51  lies  off  scale  at  a  ratio  of  2.6  x  10'5  and  an  interference  range  of  1.8  x 
104  km.  This  evaluation  as  well  as  that  for  figure  50  bears  little  resemblance  to  the 
other  evaluations.  This  result  reinforces  our  previous  conclusion  that  the  coupling  char¬ 
acteristics  for  S  =  0,  in  which  p  is  designed  to  superimpose  a  pair  of  vertical  asymp¬ 
totes,  are  radically  different  from  those  for  critical  frequencies.  In  contrast,  the  cou¬ 
pling  characteristics  for  S  =  0  and  p  =  1  for  the  free-free  and  rigid-rigid  boundary 
conditions  are  rather  similar  to  those  for  critical  frequencies. 
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AMPLITUDE  RATIO 


Figure  52.  Correlation  between  the  mode 
interference  range  and  the  ratio  of 
minimum  to  maximum  amplitude  of  the 
lower  order  mode. 


H.  EIGENVALUE  DIFFERENCES  AT  ESTIMATES  OF  CRITICAL 
FREQUENCY 

This  section  presents  numerical  results  that  follow  from  section  2J.  In  the  early 
stages  of  this  report,  we  examined  various  cases  of  p  and  S  in  a  somewhat  haphazard 
manner.  These  cases  indicated  that  the  value  of  y  at  the  critical  frequency  was  a  more 
important  parameter  than  either  the  parameters  p  or  S  from  which  y  is  derived.  Our 
initial  impression  was  that  when  y  was  greater  than  about  -1,  there  was  no  significant 
coupling  and  that  when  y  was  more  negative  than  about  -2,  the  modes  were  too 
closely  coupled  to  be  of  interest.  Thus,  it  initially  appeared  that  the  cases  of  most 
interest  fell  roughly  between  y  =  -1  and  y  =  -2.  Figure  14  is  a  good  example  of  this 
situation.  Section  2J  was  developed  as  a  analysis  tool  to  determine  if  this  impression 
was  correct. 
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Figure  53  presents  the  results  of  the  analysis  for  the  free-rigid  boundary  condition. 
The  results  only  apply  to  modes  1  and  2.  The  vertical  scale  is  the  estimate  of  the 
difference  between  eigenvalues  as  given  by  the  discriminate  of  eq.  A-6.  The  horizontal 
scale  is  p.  Each  curve  corresponds  to  a  fixed  value  of  y  ranging  from  -0.1  to  -10.  The 
procedure  for  generating  a  curve  is  to  solve  eq.  107  for  S  in  terms  of  y  and  p  for  m 
and  n  =  1.  We  take  the  fixed  value  of  y  and  evaluate  S  for  some  choice  of  p.  For  each 
set  of  y,  p,  and  S  we  evaluate  the  discriminate  using  eqs.  101  and  102.  Each  curve  is 
generated  by  stepping  through  the  desired  range  of  p. 


p 

Figure  53.  Estimate  of  eigenvalue  separation  as  a  function  of  p 
for  various  values  of  y  for  the  estimated  critical  frequency  for 
the  free-rigid  boundary. 


Figures  54  and  55  are  the  respective  counterparts  of  figure  53  for  the  free-free  and 
rigid-rigid  boundary  conditions.  Here,  we  use  eqs.  110  or  114  to  solve  for  S  and  evalu¬ 
ate  the  discriminant  using  eqs.  101  or  102. 
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We  note  that  the  parameters  S  is  not  shown  in  figures  53  to  55.  In  the  case  of 
figures  54  and  55,  p  =  1  corresponds  to  S  =  0.  Thus,  for  p  >  1,  S  is  positive  and 
increases  with  increasing  p.  For  p  <  1,  S  is  negative  and  becomes  more  negative  as  p 
goes  to  zero.  In  the  case  of  figure  53,  p  =  0.6601015  corresponds  to  S  =  0  and  S  is 
positive  for  p  >  0.6601015  and  negative  for  p  <  0.6602015. 

We  note  that  for  values  of  p  greater  than  about  0.2  or  0.3  the  discriminant  de¬ 
creases  as  y  becomes  more  negative.  This  supports  our  original  impression  that  y  is 
the  important  parameter  and  that  the  mode  separation  increases  as  y  becomes  less 
negative.  For  values  of  p  below  0.2  or  0.3  the  curves  of  y  are  no  longer  nested  but 
cross  over  each  other.  However  the  evaluations  in  figures  53  to  55  in  this  area  are 
incomplete  and  need  further  work  beyond  the  scope  of  this  report.  The  chief  purpose 
of  figures  53  to  55  is  to  serve  as  a  guide  in  selecting  the  choices  p  of  and  S  to  be 
used  as  examples  of  the  various  patterns  the  eigenvalue  solutions  take  on  as  presented 
in  the  next  section. 

I.  EIGENVALUES  FOR  VARIOUS  COMBINATIONS  OF  p  AND  S 

This  section  presents  canonical  eigenvalues  and  associated  sound-speed  profiles  for 
various  combinations  of  p  and  S.  This  is  not  a  thorough  investigation,  but  rather  a 
sampling  to  give  a  feel  for  the  eigenvalue  patterns  for  a  variety  of  choices  of  p  and  S. 

We  start  the  analysis  with  an  examination  of  small  negative  values  of  y  to  deter¬ 
mine  if  it  is  possible  to  have  close  coupling  for  critical  points  close  to  y  =  0.  Figure  56 
presents  the  eigenvalues  for  S  =  0.8948  and  p  =  3.9  for  the  free-rigid  boundary.  These 
parameters  correspond  to  the  last  point  evaluated  on  the  y  =  -0.2  curve  in  figure  53. 
The  value  of  the  discriminate  is  0.1075. 

We  see  that  the  separation  between  mode  eigenvalues  can  indeed  be  very  small  for 
critical  points  with  small  negative  values  of  y.  This  not  only  occurs  for  the  critical 
points  for  modes  1  and  2  but  also  for  modes  4  and  5  at  a  value  of  y  even  closer  to 
zero  than  y  =  -0.2.  This  separation  is  also  small  for  the  critical  points  for  modes  2 
and  3  and  3  and  4  for  which  y  is  closer  to  zero  than  y  =  -0.5. 

Figure  57  presents  the  eigenvalues  for  S  =  0.8574  and  p  =  3.8  for  the  free-free 
boundary.  These  parameters  correspond  to  the  last  point  evaluated  on  the  y  =  -  0.2 
curve  in  figure  54.  The  value  of  the  discriminate  is  0.108.  Small  differences  occur  for 
modes  1  and  2  at  y  =  -0.2  and  for  modes  2  and  3,  3  and  4,  and  4  and  5  for  y  closer 
to  zero  than  -0.5. 

Figure  58  presents  the  eigenvalues  for  S  =  0.519  and  p  =  2.1  for  the  rigid-rigid 
boundary.  In  this  case,  we  did  not  take  the  endpoint  of  the  y  =  -0.2  curve  in  figure  55 
because  the  discriminate  was  significantly  larger  than  0.1.  Here,  we  selected  p  =  2.1  as 
the  point  where  the  y  =  -0.5  curve  leads  to  a  discriminate  of  0.104.  Here  the  small 
differences  occur  for  modes  1  and  2  for  y  =  -0.5  and  for  modes  2  and  3,  3  and  4,  and 
4  and  5  for  values  of  y  between  about  -1  and  -2.5. 
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Figure  56.  Eigenvalues  for  S  =  0.895,  p  =  3.9, 
free-rigid  boundary. 


Figure  57.  Eigenvalues  for  S  =  0.857,  p  =  3.8, 
free-free  boundary. 


X 

Figure  58.  Eigenvalues  for  S  =  0.519,  p  =  2.1, 
rigid-rigid  boundary. 

Figures  59  a  to  c  present  the  sound-speed  profiles  associated  with  the  parameters 
of  figures  56  to  58  respectively.  As  usual,  the  standard  parameters  for  the  surface 
layer  were  used.  These  profiles  are  examples  approaching  the  limiting  case  of  a  posi¬ 
tive-gradient  layer  overlapping  an  isospeed  layer.  In  this  limit,  S  goes  to  1  and  p  goes 
to  infinity.  In  this  limit,  eqs.  42  or  44  reduce  to  y  =  0.  Thus,  in  the  limit  the  slant 
asymptotes  of  figures  56  to  58  collapse  to  a  horizontal  line  at  y  =  0.  This  result  sug¬ 
gests  that  close  coupling  between  modes  can  occur  for  very  small  negative  values  of  y. 
This  is  done  by  letting  S  approach  closer  to  1  and  p  approach  large  positive  values. 

Our  next  six  examples,  were  chosen  to  determine  if  relatively  large  negative  values 
of  y  always  lead  to  extremely  close  coupling  between  modes  or  if  modest  mode  cou¬ 
pling  can  occur.  The  parameters  p  and  S  for  all  six  examples  where  chosen  so  that  y  = 
-5  for  the  intersection  between  the  asymptotes  associated  with  modes  1  and  2.  For  the 
first  three  examples,  p  was  chosen  in  figures  53  to  55  so  that  the  discriminant  would 
be  at?'  i  0.1.  Figure  60  is  for  S  =  -0.549  and  p  =  0.4  for  the  free-rigid  boundary. 
Figure  61  is  for  S  =  -2.793,  p  =  0.32  for  the  free-free  boundary.  Figure  62  is  for 
S  =  -1.003,  p  =  0.38  for  the  rigid-rigid  boundary.  Only  the  A  mode  eigenvalues  are 
shown  here  as  will  be  the  case  for  other  examples  to  be  presented  later.  In  all  three 
figures,  modes  1  and  2  have  extremely  close  coupling,  which  is  not  comensurate  with 
our  discriminate  value  of  0.1.  More  will  be  said  of  this  later.  However,  a  fortuitous 
result  in  figure  61  shows  a  significant  displacement  between  modes  2  and  3  at  the 
point  of  closest  approach  for  a  value  of  y  somewhat  more  negative  that  -2. 
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Figure  60.  Eigenvalues  for  S  = 
p  =  0.4,  free-rigid  boundary. 
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Figure  62.  Eigenvalues  for  S  =  -1.003, 
p  =  0.380,  rigid-rigid  boundary. 


This  latter  example  suggests  that  significant  displacements  can  occur  for  modest 
values  of  negative  y.  For  the  next  three  examples,  p  was  chosen  from  figures  53  to  55 
to  yield  large  values  of  the  discriminate,  viz,  0.8,  0.6,  and  about  0.6.  The  idea  here 
was  to  determine  if  we  could  separate  modes  1  and  2  at  the  critical  frequency  by 
increasing  the  discriminant  of  our  approximation.  Figure  63  is  for  S  =  -2.991,  p  =  0.2 
for  the  free-rigid  boundary.  Figure  64  is  for  S  =  -29.93,  p  =  0.1,  for  the  free-free 
boundary.  Figure  65  is  for  S  =  -4.062,  p  =  0.2  for  the  rigid-rigid  boundary.  Figure  63 
displays  a  significant  displacement  between  modes  1  and  2  for  y  between  -3  and  -4. 
Figure  64  displays  extremely  close  coupling  between  modes  1  and  2.  However,  there  is 
a  relatively  large  displacement  between  modes  3  and  4  for  y  between  -3  and  -4.  Fig¬ 
ure  65  displays  a  small  but  definite  displacement  between  modes  1  and  2  for  y  be¬ 
tween  -3  and  -4. 

Figures  63  to  65  clearly  indicate  why  the  approximation  for  the  displacement  be 
tween  modes  at  the  critical  frequencies  fail.  The  main  problem  is  that  the  estimates  of 
critical  frequencies  of  section  2F  given  by  the  intersections  of  the  asymptotes  are  not 
good.  The  critical  frequencies  of  figures  63  ,  65,  and  64  lie  well  above  the  intersection 
at  y  =  -5.  This  is  also  true  to  a  somewhat  lesser  extent  for  figures  60  to  63.  Indeed, 
there  are  some  good  indications  that  the  expansions  of  Section  2J  are  reasonably  good 
if  we  examine  the  mode  displacements  at  y  =  -5  rather  than  at  the  critical  frequency. 
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Figure  63.  Eigenvalues  for  S  =  -2.991,  p  =  0.2, 
free-rigid  boundary. 
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Figure  64.  Eigenvalues  for  S  =  -29.93,  p  =  0.1, 
free-free  boundary. 
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Figure  65.  Eigenvalues  for  S  =  -4.062,  p  =  0.2, 
rigid-rigid  boundary. 

For  example,  in  figure  63  this  displacement  is  about  0.8  as  was  originally  determined 
from  figure  53.  Similarly  the  mode  displacements  at  y  =  -5  are  large  for  figures  64 
and  65  and  significantly  smaller  for  figures  60  to  62  in  line  with  our  original  choice  of 
displacement. 

The  sound-speed  profiles  associated  with  figures  60  to  65  are  presented  in  figures 
66  and  67.  These  are  the  first  examples  of  negative  values  of  S.  Negative  S  means  that 
the  bottom  boundary  sound  speed  is  smaller  than  the  surface  sound  speed.  The  more 
negative  the  value  of  S,  the  lower  the  bottom  sound  speed.  Figure  67b  for  S  =  -29.930 
is  the  most  extreme  case  shown  here  with  bottom  sound  speeds  some  50  odd  m/s 
below  the  surface  sound  speed.  These  negative  values  of  S  give  positive  slopes  on  the 
slant  asymptotes  of  figures  60  to  65  as  given  by  eqs.  42  or  43.  As  the  bottom  sound 
speed  goes  to  zero  as  a  limit  S  goes  to  -  and  the  slope  of  the  slant  asymptotes  goes 
to  1  as  a  limit.  This  limit  is  probably  related  to  the  unity  slope  for  the  line  y  =  x. 
Asymptotes  with  slopes  between  -h»  and  one  intersect  the  line  y  =  x  for  positive  values 
of  y.  These  have  no  significance  as  this  is  beyond  the  regime  of  mode  coupling.  How¬ 
ever,  suppose  the  slope  of  the  slant  asymptotes  could  be  less  than  one.  These  slant 
asymptotes  would  intersect  the  line  y  =  x  for  large  negative  values  of  y.  This  intersec¬ 
tion  would  represent  an  impossible  physical  situation  because  y  =  x  represents  zero 
frequency,  whereas  large  negative  y  represents  the  large  end  of  the  frequency  spec¬ 
trum. 
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We  note  another  consequence  of  negative  S.  At  sufficiently  large  frequencies,  the 
low-order  modes  will  be  associated  with  the  slant  asymptotes.  (This  is  just  the  opposite 
of  the  case  of  positive  S  where  at  high  frequencies,  the  low-order  modes  are  associated 
with  the  vertical  asymptotes).  The  reason  is  that  for  negative  S  the  modes  are  initially 
trapped  in  the  lower  duct,  because  the  sound  speed  at  the  bottom  boundary  is  smaller 
than  that  at  the  surface. 

For  our  next  six  examples,  we  chose  a  more  modest  value  of  y  =  -2  with  values  of 
S  slightly  positive  in  figures  68  to  70  and  slightly  negative  in  figures  71  to  73.  In  fig¬ 
ure  68  the  value  of  p  =  0.73  was  chosen  from  figure  53  for  a  discriminant  of  0.1.  In 
figure  69  p  =  1.05  was  chosen  from  figure  54  for  a  discriminant  of  0.015.  In  figure  70 
p  =  1.05  was  chosen  from  figure  55  for  a  discriminant  of  0.025.  The  values  of  S  for 
these  three  cases  is  positive.  We  see  that  modes  1  and  2  are  strongly  coupled  in  fig¬ 
ures  69  and  70  corresponding  to  small  values  of  the  discriminant,  whereas  figure  68 
displays  a  separation  comensurate  with  a  larger  discriminant  of  0.1.  In  figure  71  the 
value  of  p  =  0.57  was  chosen  from  figure  53  for  a  discriminant  of  0.2.  In  figure  72 
p  =  0.56  was  chosen  from  figure  54  for  a  discriminant  of  0.1.  In  figure  73  p  =  0.73 
was  chosen  from  figure  55  for  a  discriminant  of  0.1.  The  values  of  S  for  these  3  cases 
is  negative.  Here,  we  see  the  most  separation  between  modes  1  and  2  in  figure  71  with 
about  one-half  this  amount  of  separation  in  figures  72  and  73.  Figures  72  and  73 
contain  fortuitous  results  in  the  relatively  large  separation  between  modes  4  and  5  for 

y  =  -1. 

The  sound-speed  profiles  associated  with  figures  68  to  73  are  presented  in  figures 
74  and  75.  These  profiles  are  characterized  by  surface  and  bottom  sound  speeds  rea¬ 
sonably  close  to  each  other.  They  are  not  nearly  as  extreme  as  the  profiles  of  figures 
66  and  67. 

In  closing  this  section,  we  believe  that  the  results  of  figures  71  to  73  answer  our 
original  speculation  about  large  separations  between  modes  for  relatively  large  negative 
values  of  y.  The  key  is  to  pick  small  negative  values  of  S  together  with  values  of  p 
leading  to  the  desired  y.  It  would  be  interesting  to  see  if  this  approach  will  lead  to 
large  separations  for  even  more  negative  values  of  y.  However,  this  would  entail  more 
extensive  and  careful  evaluation  in  figures  53  to  55,  which  is  beyond  the  scope  of  the 
present  report. 
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Figure  68.  Eigenvalues  for  S  =  0.0983,  p  =  0.73, 
free-rigid  boundary. 


Figure  69.  Eigenvalues  for  S  =  0.0501,  p  =  1.05, 
free-free  boundary. 
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Figure  70.  Eigenvalues  for  S  =  0.0314,  p  =  1.05, 
rigid-rigid  boundary. 
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Figure  71.  Eigenvalues  for  S  =  -0.184,  p  =  0.57, 
free-rigid  boundary. 
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Figure  72.  Eigenvalues  for  S  =  -1.18,  p  =  0.56, 
free-free  boundary. 
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Figure  75.  Sound-speed  profiles  associated  with 
figures  71  to  73  respectively. 
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4.  AREAS  FOR  FURTHER  INVESTIGATION 


This  brief  section  points  up  some  of  the  areas  requiring  additional  investigation. 

One  of  the  primary  areas  of  interest  is  the  treatment  of  the  normalization  process. 
Although  the  present  approach  is  correct,  it  is  not  satisfactory  to  have  to  deal  with  Dn 
that  ranges  from  small  to  huge  or  even  divergent  values.  An  investigation  of  alternative 
methods  of  defining  the  unnormalized  eigenfunctions  so  that  Dn  is  well  behaved  is 
recommended.  Some  of  these  alternative  methods  are  suggested  in  section  3A. 

Another  area  of  interest  is  to  determine  the  common  curve  Dca  approaches  for  all 
modes  at  high  frequencies  in  figure  31  and  also  with  some  additional  complication  in 
figure  20.  A  related  investigation  is  why  Dca,  for  the  rigid-rigid  boundary  condition  in 
figure  38,  does  not  approach  a  common  curve  for  large  frequencies  as  is  the  case  for 
other  boundary  conditions. 

Some  aspects  of  the  example  of  section  3F  still  pose  an  interesting  enigma  in  that 
some  features  associated  with  coupled  modes  are  present  such  as  almost  identical 
amplitudes  in  the  lower  duct  but  radically  different  amplitudes  in  the  upper  duct. 
Although  a  detailed  investigation  of  the  mathematics  of  this  example  indicate  that  the 
results  are  correct,  we  have  no  rationale  from  a  physical  perspective  as  to  why  this 
differs  from  other  coupled  mode  configurations.  It  would  also  be  of  interest  to  deter¬ 
mine  if  there  is  some  other  approach  to  eigenvalues  that  would  not  lead  to  the  diver¬ 
gent  values  of  Dn  encountered  for  this  special  case. 

It  would  be  interesting  to  make  a  thorough  investigation  of  the  accuracy  of  the 
mode  displacement  at  the  critical  frequencies  as  given  by  the  approximations  of  section 
J.  The  present  work  suggests  that  accuracy  may  be  quite  good  for  positive  values  of  S 
and  slightly  negative  values  of  S.  The  chief  problem  is  with  large  negative  values  of  S, 
where  the  critical  frequencies  do  not  lie  close  to  the  approximation  given  by  the  inter¬ 
section  of  asymptotes.  It  would  be  desirable  to  formulate  a  new  approach  for  predict¬ 
ing  critical  frequencies  for  large  negative  values  of  S,  although  it  is  not  obvious  how  to 
proceed  in  the  formulation. 

It  would  also  be  desirable  to  clean  up  the  curves  of  figures  53  to  55  for  large  nega¬ 
tive  values  of  y.  This  would  be  a  first  step  in  evaluating  our  premise  that  large  dis¬ 
placements  may  be  found  for  large  negative  y  with  the  use  of  small  negative  values  of 
S  and  appropriate  values  of  p. 

In  the  early  stages  of  this  report,  we  spent  some  time  investigating  the  mathemati¬ 
cal  principle  involved  in  section  3E  in  which  the  eigenvalues  for  one  mode  remains 
almost  stationary,  whereas  those  of  the  other  mode  move  rapidly  as  p  is  moved  above 
or  below  the  values  leading  to  common  vertical  asymptotes.  Here,  the  physical 
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principle  is  related  to  an  abhorrence  of  degenerate  eigenvalues.  However,  the  detailed 
mathematical  whys  and  wherefores  remain  of  interest.  We  did  not  complete  this  inves¬ 
tigation.  However,  the  preliminary  results  appeared  promising  and  should  be  pursued 
further.  Our  approach  was  to  investigate  the  quadratic  expression  for  coincident  roots. 
Here,  the  discriminant  goes  to  zero  leaving  only  one  value  no  matter  which  sign  is 
picked  for  the  square  root  of  the  discriminant.  However,  when  we  differentiate  the 
expression  with  respect  to  the  parameter  p,  the  derivative  was  different  for  the  choice 
of  signs  on  the  square  root  of  the  discriminant.  Here,  one  choice  of  signs  leads  to  a 
rate  of  change  with  respect  to  p  of  zero  whereas  the  other  choice  leads  to  a  rate  of 
change  equal  to  twice  that  of  a  simple  zero. 

In  closing  this  report,  we  note  that  we  have  achieved  our  original  goal  in  formulat¬ 
ing  the  canonical  eigenvalue  approach.  This  goal  was  to  reduce  the  number  of  parame¬ 
ters  (in  this  case  from  six  to  two)  so  that  the  salient  features  of  double-duct  propaga¬ 
tion  could  be  analyzed.  However,  even  this  simplification  leaves  a  somewhat 
formidable  problem  in  that  the  investigation  still  involves  the  behavior  of  eigenvalues 
in  the  two-dimensional  space  of  S  and  p.  We  attempted  to  solve  this  problem  by  intro¬ 
ducing  the  parameter  y  for  the  intersection  of  the  asymptotes  for  S,  p,  and  the  bound¬ 
ary  conditions.  We  found  that  y  does  indeed  appear  to  provide  more  information  about 
the  nature  of  the  eigenvalues  than  does  p  or  S  alone.  The  introduction  of  the  dis¬ 
criminant  in  the  representation  of  figures  53  to  55  also  helped  in  the  analysis.  How¬ 
ever,  these  approaches  are  not  adequate,  particularly  for  large  negative  values  of  S. 

We  have  come  a  long  way,  but  there  is  a  long  way  to  go  before  a  satisfactory  under¬ 
standing  is  available  for  this  most  simple  of  double  ducts  to  say  nothing  of  more  com¬ 
plicated  double-duct  configurations. 
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APPENDIX  A:  TAYLOR  SERIES  EXPRESSIONS  ABOUT  THE 
ASYMPTOTE  INTERSECTIONS 


In  this  appendix,  we  estimate  the  horizontal  separation  between  mode  eigenvalues 
at  the  point  where  the  asymptotes  intersect.  This  separation  is  of  interest  because  it 
indicates  the  amount  of  coupling  between  the  modes  at  the  critical  frequencies. 

Let  the  eigenvalue  equation  be  given  as 


G  (x,  y)  =  0  .  (A-l) 

We  expand  G  in  a  Taylor  series  about  the  intersection  point  (x0,  y0).  We  shall  expand 
eq.  A-l  for  the  case  of  y  =  y0  so  as  to  deal  with  a  one-dimensional  rather  than  a 
two-dimensional  expansion.  Equation  A-l  is  also  a  function  of  S,  p,  mode  numbers, 
and  boundary  conditions,  all  of  which  are  considered  fixed  and  enter  into  the 
expansions  as  parameters.  The  first  three  terms  of  the  expansion  of  eq.  A-l  may  be 
written  as 

G(x,  y0)  =  G(x0,  yG)  +  (x  -  xQ)  <5G/<5x  +  (x  -  x0)22_1<52G/<5x2  .  (A-2) 

The  roots  of  this  quadratic  expression  set  to  zero  may  be  written  as 


(x  -  x0)  =  [-  b  ±  (b2  -  4c)1/2]/2  (A-3) 

where 

b  =  2  (<5G/<5x)  /  (<32G/(5x2)  (A-4) 

and 

c  =  2G(x0,  y0)/(<52G/<5x2)  .  (A-5) 

The  displacement  between  the  roots  is 

X2  -  xj  =  (b2  -  4c)1/2  .  (A-6) 

Equation  A-6  then  is  an  approximation  to  the  horizontal  separation  between  the 
eigenvalues. 

The  next  step  is  to  derive  expressions  for  <52G/(Sx2-  We  differentiate  eq.  29  for  the 
free-rigid  boundary  with  respect  to  x  and  with  the  use  of  Eqs.  35  nd  32  obtain 


A-l 


(A-7) 


62G\s/dx2  =  -  [x  +  £2oP4(l  -  S)2]Gi5 

-2£20p2(1-S)G18  +  p4(1-S)2G17. 

To  obtain  the  counterparts  of  eq.  A-7  for  other  boundary  conditions,  we  need  the 
expression 


<5G18/(5x  =  xGn  -  p2(l  -  S)G16  .  (A-8) 

We  differentiate  eq.  32  for  the  free-free  boundary  condition  with  respect  to  x  and 
with  the  use  of  eqs.  A-8  and  29  obtain 

d2G  17/<Sx2  =  -  [x  +  £20p4(1  -  S)2]G17  +  2p2(l  -  S)G16  .  (A-9) 

We  differentiate  eq.  35  for  the  rigid-rigid  boundary  condition  with  respect  to  x  and 
with  the  use  of  eqs.  29  and  A-8  obtain 

(52G16/(5x2  =  G15  -  [x  +  Up4(1  -  S)2]G16  +  2x£20p2(1  -  S)  G17  +  xp4(1  -  S)2G18  .  (A-10) 

We  now  evaluate  the  dominant  term  of  the  factors  entering  into  the  expansions  of 
the  eigenvalue  equation.  The  largest  term  for  large  negative  y  is  G4.  The  smallest  term 
is  G,.  The  terms  G2  and  G3  fall  between  and  are  of  equal  importance. 

In  the  case  of  the  free-rigid  boundary  condition  Ai(-xo)  =  0  and  Ai'(-S2o)  =  0.  The 
dominant  terms  then  become 


Gi5(x0,  y<>)  =  -Bi(-x)  Bi'(-  £20)Gi  (A-ll) 

(5Gls/<5x  =  -  Ai'(-x)Bi'(-  £20)G2  +  C2oP2(1  -  S)Bi(-  x)Ai(-  £20)G3  (A-12) 

<52Gls/<5x2  =  2^20p2(1  -  S)  Ai'(-  x)  Ai(-  £20)  G4  .  (A-13) 

The  counterparts  of  eqs.  A-ll  to  A-13  for  the  free-free  boundary  conditions  are: 
G17(x0,  y0)  =  -  Bi(-  x)  Bi(-  U)  G,  (A-14) 

6Gv/6x  =  -  Ai'(-  x)Bi(-  £2o)G2  -  P2(l  -  S)Bi(-  x)Ai'(-  ^20)G3  (A-15) 

and 

62G\i/dx2  =  -  2p2(1  -  S)  Ai'(-  x)  Ai'(-  C20)G4  .  (A-16) 


A-2 


The  counterparts  of  eqs.  A-ll  to  A-13  for  the  rigid-rigid  boundary  conditions  are 


G16(x0,y0)  =  -Bi'(-x)  Bi’(-ko)G,  (A-17) 

c5G16/(Sx  =  xAi(-  x)Bi'(-  ?2o)G2  +  ?20P2(1  -  S)Bi’(-  x)Ai(-  J2 o)G3  (A-18) 

and 

<5G16/<5x2  =  -  2x£20p2(1  -  S)  Ai(-  x)  Ai(-  £20)G4  .  (A-19) 

Next,  consider  the  evaluation  of  eqs.  A-4  and  A-5  for  the  free-rigid  boundary 
condition.  The  substitution  of  eqs.  A-ll  to  A-13  leads  to 

b  =  -  Bi'(-  £2o)(G2/G4)/koP2(l  -  S)  Ai(-  £20)  +  Bi(-  x)(G3/G4)Ai'(-  x)  (A-20) 

and 

c  =  -  Bi(-  x)Bi'(-  £20) (G i / G4) / ^20 p2 ( 1  -  S)Ai'(-  x)Ai(-  £20)  .  (A-21) 

We  may  eliminate  the  Bi  and  Bi'  functions  in  eqs.  A-20  and  A-21  by  means  of 

Bi(- x)  =  -  [tt  Ai'(- x)]"1  (A-22) 

and 

Bi'(-  £20)  =  —  [tt  Ai(-  ^o)]'1  (A-23) 

Equation  A-22  and  A-23  follow  from  the  Wronskian  evaluated  at  various  zeros  of  the 
Airy  functions.  Equations  101  to  105  of  the  main  text  then  follow  from  eq.  A-20  to 
A-33. 

For  the  free-free  boundary  condition,  eqs.  A-14  to  A-16  lead  to 
b  =  -  Bi'(-  ?2o)(G2/G4)/p2(1  -  S)  Ai'(-  C20)  +  Bi(-  x)(G3/G4)Ai'(-  x)  (A-24) 

and 

c  =  -  Bi(-  x)Bi(-  ?2o)(G1/G4)/p2(1  -  S)Ai'(-  x)Ai(-  £20)  .  (A-25) 

Equations  118  and  109  of  the  main  text  then  follow  from  eqs.  A-22  to  A-25. 

For  the  rigid-rigid  boundary  condition  eq.  A-17  to  A-19  lead  to 

b  =  -  Bi'(-  S2o)(G2/G4)/£20p2(1  -  S)  Ai(-  £20)  -  Bi'(-  x)(G3/G4)/xAi'(-  x)  (A-26) 
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and 


c  =  Bi'(-  x)  Bi'(-  C20) (G1/G4)/x^20p2(1  -  S)  Ai'(-  x)  Ai(-  £20)  .  (A-27) 

Equations  111  and  112  of  the  main  text  then  follow  from  eqs.  A-26,  A-27,  A-22,  and 
A-23. 
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APPENDIX  B:  EXTENSION  OF  THE  CANONICAL 
EIGENVALUE  METHOD  TO  INCLUDE  DISCONTINUITIES  OF 
SOUND  SPEED  AND  TO  INCLUDE  THE  EFFECTS  OF 
LAYERS  WITH  DIFFERENT  CONSTANT  DENSITIES 


Reference  B1  presents  some  recent  work  by  Computer  Sciences  Corporation  (CSC) 
under  contract  to  NOSC,  Code  541.  Portions  of  this  work  deal  with  updating  the 
canonical  eigenvalue  method  so  as  to  be  able  to  treat  models  of  bottom  interaction 
which  involve  sound-speed  discontinuities,  layers  with  different  densities,  and  the  use 
of  complex  parameters  to  model  attenuation  in  bottom  sediments.  Since  there  is  no 
immediate  plan  to  publish  reference  Bl,  it  seemed  advisable  to  include  the  new  work 
on  the  canonical  eigenvalue  method  here  in  the  form  of  this  appendix.  Some  of  the 
models  of  bottom  interaction  involve  double-duct  propagation.  However,  the  numerical 
evaluation  of  these  cases  will  not  be  carried  out. 


1.  MODIFICATIONS  TO  THE  CANONICAL  EIGENVALUE  METHOD  FOR 
SOUND-SPEED  DISCONTINUITIES 

To  use  the  canonical  eigenvalue  method  in  the  preparation  of  controls  for  subbot¬ 
tom  propagation,  it  is  necessary  to  modify  the  treatment  of  reference  1  to  include  dis¬ 
continuities  in  the  sound  speed  at  layer  interfaces.  Let  the  sound  speed  at  the  bottom 
of  layer  i  be  given  by  Q+io  and  at  the  top  of  layer  i+1  be  given  by  Ci+li].  When  the 
sound  speed  is  continuous  at  the  interface  between  layer  i  and  i+1, 

Ci+1,0  =  Cj+1,1  • 

Our  interest  here  is  in  the  case  where  eq.  B-l  does  not  apply.  Each  interface  with  a 
sound-speed  discontinuity  introduces  one  more  parameter  into  the  profile  representa¬ 
tion. 

From  eq.  4  of  reference  1,  we  determine  that 

-a,/ai+i  =  (-yi0/Vi+i)'/3  Ci+M/Ci+1,o  =  Pi  .  (B-2) 

The  first  of  our  new  canonical  profile  parameters  is  P;  defined  by  eq.  B-2.  We  note 
that  when  eq.  B-l  holds, 


Pi 


Pi  • 


(B-3) 
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Our  next  step  is  to  determine  Mi+1  and  Ni+J  in  the  counterpart  of  eq.  18  of  refer¬ 
ence  1.  This  counterpart  is 

£1+1,0  =  Mi+1  £ij  +  Ni+i  £i0  .  (B-4) 

We  proceed  in  a  similar  manner  to  the  development  of  eqs.  18  to  31  of  reference  1. 
The  counterparts  of  eqs.  30  and  31  become 

Ni+1  =  Pi2Si+1  (B-5) 

and 

Mi+1=  P.MI-Sm)  (B-6) 

where 

s,„  =  [(Cil/C„2,o)2  -  l]/l(C„/Ci„.o)2  -  1]  (B-7) 

Equation  B-7  is  the  counterpart  of  eq.  29  of  reference  1  and  represents  the  second  of 
our  new  canonical  profile  parameters.  When  eq.  B-l  applies,  Ni+1,  Mi+1,  and  Si+1 
reduce  respectively  to  Ni+i,  Mi+1,  and  Si+i,  of  reference  1.  We  note  that  eq.  B-7 
involves  the  sound  speeds  at  the  top  and  bottom  of  layer  i  and  at  the  bottom  of 
layer  i+1. 

Next,  we  must  develop  an  expression  for  a  third  canonical  profile  parameter,  which 
is  necessary  to  account  for  the  extra  third  profile  parameter  for  a  layer  with  sound- 
speed  discontinuity.  We  assume  that  the  Airy  function  argument  at  the  top  of  layer  i+1 
may  be  written  an  a  linear  combination  of  the  arguments  at  the  top  and  bottom  of  lay¬ 
er  i.  Thus, 

£i+i ,  1  -  Mi+1  +Ni+j  £i0  .  (B-8) 

Our  task  is  to  determine  Mi+1  and  Ni+1.  We  proceed  in  a  manner  similar  to  that  of 
eqs.  18  to  31  of  reference  1.  It  is  straightforward  to  determine  expressions  for  Mi+1 
and  Ni+1.  However,  an  appropriate  third  canonical  parameter  was  not  all  apparent  from 
these  original  expressions.  By  analogy  with  eqs.  B-5  and  B-6  we  assumed  that  the  third 
parameter  was  Si+1  and  that 


We  then  used  our  original  expressions  for  Ni+1  or  Mi+1  to  solve  for  Si+1  to  obtain 

Sw  =  KCi^.o/Cui.i^KCuM/Qti^-lJ/fQ^.o/Q,)2-!]  .  (B-ll) 

This  then  is  the  third  canonical  profile  parameter.  When  eq.  B-l  applies,  Si+1  =  1, 

Mi+1  =  0,  and  Ni+1  =  i2.  Thus,  eq.  B-8  reduces  to 

£i+l,l  =  P  i  2  £i0  . 

which  is  the  result  of  eq.  10  of  reference  1  for  the  continuous  profile.  We  note  that 
eq.  B-ll  involves  the  sound  speeds  at  the  top  and  bottom  of  layer  i  and  at  the  top  of 
layer  i+1. 

We  note  that  the  original  expressions  for  Ni+,  or  Mi+1  contained  Pi 2,  but  not  Pi  2. 
The  (Ci+i.o/Ci+i,])2  factor  was  inserted  in  eq.  B-ll  so  that  Pj2  rather  than  P2 
appears  in  eqs.  B-9  and  B-10.  This  allowed  us  to  define  Si+1  such  that  Ni+1  and  Mi+1 
can  be  expressed  as  simple  functions  of  the  canonical  parameters. 

Equations  B-4  and  B-8  or  their  continuous  counterparts  allow  us  to  determine  all 
Airy  function  arguments  as  a  linear  function  of  x  and  y  by  a  recursive  process  similar 
to  that  discussed  in  eqs.  32  to  39  of  reference  1. 

Consider  a  profile  consisting  of  m  boundaries,  ni  interfaces  without  discontinuity 
and  n2  interfaces  with  sound  speed  discontinuity.  The  total  number  of  profile  parame¬ 
ters  is  2n,  +  3n2  +  m  +  1  parameters.  Here,  we  specify  n,  values  of  Pi,  n,  +  m  -2  val¬ 
ues  of  Si+i,  n2  values  of  Pi  2  and  n2  values  each  of  Si+1  and  Si+]  for  a  total  of 
2n,  +  3n2  +  m  -2  canonical  parameters.  If  we  then  select  C1f  C20,  and  Yt  we  have  the 
requisite  number  of  parameters.  These  three  profile  parameters  plus  the  set  of  canoni¬ 
cal  parameters  can  be  used  to  generate  the  sound-speed  profile. 

Assume  that  the  first  interface  with  a  discontinuous  sound  speed  is  at  the  bottom  of 
layer  i  and  that  all  parameters  for  layer  1  to  the  bottom  of  layer  i  have  been  generated 
by  the  recursive  method  of  reference  1.  The  first  step  is  to  use  the  given  value  of  Si+1, 
the  known  values  of  and  Cn  and  Cit)io  to  generate  Ci+]ji  from  eq.  B-ll.  The  second 
step  is  to  use  the  given  values  of  Pj  and  the  known  values  of  Yl0,  Ci+]>i,  and  Ci+ii0  to 
generate  "yui  from  eq.  B-2.  The  third  step  is  to  use  the  given  value  of  Si+1  and  the 
known  values  of  Cn  and  Ci+i,0  to  generate  Ci+2,o  from  eq.  B-7.  The  fourth  step  is  to 
evaluate 

y>+i,o  =  ^i+i  (Ci+2,o/Ci+1..)3-  (B-l 2) 

Wc  now  know  all  the  requisite  profile  parameters  down  to  the  bottom  of  layer  i+1.  The 
four-step  process  is  then  repeated  for  each  of  the  remaining  profile  layers  to  define  the 
entire  profile. 
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Thus  far,  our  concern  has  been  with  the  generation  of  the  Airy  functions  arguments 
and  the  characterization  of  the  sound-speed  profile  in  terms  of  canonical  parameters. 
Next,  we  examine  how  the  eigenvalue  equation  and  normalization  process  are  affected 
by  the  sound-speed  discontinuities.  This  is  relatively  simple.  In  most  cases  we  only 
have  to  replace  Pi  by  P*.  However,  we  do  not  replace  Pi2  by  P2  in  the  Airy  function 
arguments  because  we  have  already  noted  how  these  arguments  change.  For  example, 
the  eigenvalue  equation  s  modified  by  replacing  Pj  in  eq.  22  of  reference  2  with  P;  at 
any  interface  where  there  are  sound-speed  discontinuities. 

Before  we  examine  the  normalization  process  of  reference  1,  we  correct  several 
critical  typographical  errors  in  reference  1.  The  most  important  of  these  corrections  is 
to  replace  af1  by  af1  on  the  right  side  of  eq.  92  and  on  the  left  side  of  eq.  96.  In 
addition,  the  6  in  eq.  33  should  be  deleted. 

A  check  of  the  normalization  process  of  reference  1  demonstrates  that  we  may 
replace  Pj  by  Pj  in  eqs.  78  to  80,  and  91.  We  may  also  replace  Pj  by  Pj  in  eq.  93 
and  Pj2by  Pj2in  eq.  96.  However,  since  eq.  10  no  longer  applies  we  must  replace  eq. 
97  by 

Dni  =  On  ((1  +  P  ?)Fi’2(-  ?i0)  t  +  P  &„)Fi2(-  £,„)]  .  (B-13) 

Equation  B-13  differs  from  the  Dni  of  eq.  97  of  reference  1  in  that  eq.  B-8  now  holds 
and  Si+i  cannot  be  expressed  as  a  simple  multiple  of  £i0.  Similarly,  in  eqs.  79  and  90 
of  reference  1  the  Airy  function  argument  of  -P;  2  £io  must  be  replaced  by  -£j+i,i. 

2.  MODIFICATIONS  TO  THE  CANONICAL  EIGENVALUE  METHOD  FOR 
LAYERS  WITH  DIFFERENT  CONSTANT  DENSITIES 

It  is  necessary  to  modify  the  treatment  of  reference  1  and  the  extensions  of  sec¬ 
tion  1  of  this  appendix  to  include  the  effect  of  a  variable  density.  Let  the  density  in 
layer  i  be  given  by  P;.  As  discussed  in  reference  Bl,  the  treatment  of  variable  density 
consists  of  three  relatively  simple  statements.  First,  it  is  necessary  that  P  U  be  con¬ 
tinuous.  Second,  the  normalization  integrand  of  eq.  83  of  reference  1  must  include  P 
as  a  multiplicative  factor.  Third,  the  propagation  loss  expression,  i.e.,  eq.  2  of  refer¬ 
ence  B2,  contains  an  additional  multiplicative  factor  of  Ps  Ph,  where  Ps  and  Ph  is  the 
density  at  the  source  and  receiver  respectively. 

The  treatment  is  facilitated  by  introducing  the  canonical  parameter  t*\  defined  by 

H  i  =  Pi/Pi+i-  (B-14) 

The  introduction  of  density  into  the  problem  has  no  effect  at  all  on  the  arguments  of 
the  Airy  functions.  Thus,  the  arguments  remain  those  values  given  in  reference  1  with 
the  modifications  for  sound-speed  discontinuities  as  discussed  in  section  1. 


B-4 


The  eigenvalue  formulation  is  modified  by  multiplying  the  top  line  of  eq.  18  of  ref¬ 
erence  2  by  M\.  Similarly,  the  right  side  of  eq.  77  of  reference  1  is  multiplied  by  Mi. 

We  note  the  continuity  of  dU/dZ  remains  unchanged  as  a  boundary  condition.  Thus, 
eq.  22'  of  reference  2  and  eq.  78  of  reference  1  remain  unchanged.  As  a  consequence, 
the  first  product  within  the  brackets  of  eqs.  79  and  80  of  reference  1  are  multiplied  by 
Mi. 


The  treatment  of  the  normalization,  Dn,  is  somewhat  more  complicated.  Equation 
89  of  reference  1  is  replaced  by 

Dn  =  a,  'ViDc  •  (B-15) 

Equation  93  of  reference  1  is  replaced  by 

Dpi  =  V  -p  j/ij  -1;  Dp]  =  1  .  (B-16) 

j  =  i 

Equations  94,  95,  98,  and  99  of  reference  1  apply  without  change  with  the  understand¬ 
ing  that  Dpi  is  evaluated  from  eq.  B-16.  The  replacement  for  eq.  B-13  or  its  refer¬ 
ence  2  counterpart,  eq.  97,  is  complicated  by  the  change  in  the  boundary  condition  in 
eq.  77  of  reference  1.  The  end  result  is  given  by 

Dni  =  Dpl[(l  +/if'  pi  3)Fi'2(-  £j0)  +  (£i0  +/*p-iCi»i..)Fi2(-  £]0)1  .  (B-17) 

We  note  that  at  interfaces  with  no  sound-speed  discontinuity  eqs.  B-16  and  B-17  apply 
with  P;  replaced  by  P; . 

3.  COMMENTS  ON  COMPLEX  EIGENVALUES 

Pages  10  and  11  of  reference  1  discuss  some  properties  of  the  complex  component 
of  the  Airy  function  arguments,  which  are  of  further  interest.  With  more  exposure  to 
the  profile  layers  with  complex  arguments,  we  conclude  that  the  statements  in  the 
paragraph  commencing  near  the  bottom  of  page  10  are  incorrect.  In  the  case  of  a  sin¬ 
gle-layer,  positive-gradient  half  space  with  complex  parameters  and  bounded  above  by 
a  free  or  rigid  surface  the  eigenvalues  are  real.  However,  this  is  not  the  case  for  multi¬ 
layer  profiles  even  with  only  one  layer  with  complex  parameters.  If  we  define  x  and  y 
as  the  top  and  bottom  values  for  the  layer  with  complex  parameters,  the  major  contri¬ 
bution  to  Im  Cp  may  indeed  be  introduced  by  the  complex  parameters  in  eq.  17  of  ref¬ 
erence  1.  However,  complex  profile  parameters  in  one  layer  lead  to  complex  canonical 
profile  parameters  which  will,  in  turn,  lead  to  complex  eigenvalue  solutions.  Thus,  the 
attenuation  associated  with  unbounded  profiles  cannot  be  completely  separated  from 
that  associated  with  complex  profile  parameters  as  was  stated  in  reference  1. 
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Equations  59  to  62  of  reference  1  present  an  analysis  of  the  imaginary  components 
of  the  Airy  function  arguments  for  real  profile  parameters.  Next,  we  extend  this  analy¬ 
sis  to  profiles  with  discontinuities  in  the  sound  speed. 

Consider  that  x  and  y  are  associated  with  a  layer  with  real  parameters.  Then,  from 
eq.  14  of  reference  1  it  follows  that  the  requirement  for  a  real  frequency  leads  to 

Im  x  =  Im  y.  (B-18) 

Indeed,  if  the  parameters  involved  in  eqs.  7  and  8  of  reference  1  are  real,  then  Cp  is 
the  only  complex  term  on  the  right  side  of  eqs.  7  and  8  and  it  follows  that 


Im  $n  =  Im£i0  .  (B-19) 

Indeed,  eq.  B-19  holds  for  any  layer  with  real  parameters.  This  means  that  the  imagi¬ 
nary  part  of  the  Airy  function  argument  is  the  same  at  the  upper  and  lower  nterfaces 
of  each  layer. 

From  eqs.  B-4  to  B-6  and  eq.  B-19,  it  follows  that 

Im  $i+li0  =  (Mi+1  +  Ni+1)  Im  ti0  =  pj  2  Im  $i0  -  (B-20) 

Similarly,  from  eq.  97  to  99  and  eq.  B-19,  it  follows  that 

Im£i+1,i  =  (Mi+1  +  Ni+1)  Im^0=Pi2ImCio  ■  (B-21) 

Thus,  as  was  the  case  for  the  continuous  profile  the  imaginary  component  of  the  Airy 
function  arguments  only  involves  Im  x  and  the  set  of  P;2. 

Assume  now  that  the  parameters  in  layer  1  are  real,  but  are  complex  in  layer  2. 
Equations  93  and  97  may  be  written  respectively  as 

£20  =  M2  Re  x  +  N2  Re  y  +  i  p}  2  Im  y  .  (B-22) 


and 


£21  =  M2  Re  x  +  N2  Re  y  +  i  Pj  2Im  y  .  (B-23) 

In  general,  ImM2  *  ImM2  and  ImN2  *  ImN2.  Thus,  Im£20  *  Im£2i  and  the  complex 
component  of  the  Airy  function  arguments  will  differ  at  the  top  and  bottom  of  the  lay¬ 
er.  However,  for  the  case  of  real  sound  speeds  but  complex  gradients,  M2,  N2,  M2, 
and  N2  are  all  real  and  Im  £20  =  Im  £2i- 
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4.  SPECIFIC  CANONICAL  EXPRESSIONS  FOR  TWO-LAYER  MODELS 


This  section  gives  the  explicit  expressions  that  are  required  to  be  evaluated  for  the 
canonical  normal  mode  solution  for  two-layer  models  of  bottom  propagation.  Figure 
B-l  is  a  schematic  indicating  the  eight-profile  parameters  that  define  the  environments. 
These  are  the  four  sound  speeds,  two  gradients,  and  two  constant  values  of  density. 
The  top  boundary  is  a  free  surface  and  the  lower  boundary  is  a  rigid  surface.  The 
parameter  C21,  C30,  and  V2  may  be  complex  so  as  to  model  a  lossy  bottom  layer. 
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Figure  B-l.  Schematic  of  the  generic  two-layer  profile  model  with 
a  sound-speed  discontinuity  and  piecewise  constant  densities. 


The  canonical  parameters  for  this  profile  class  are  given  by 


p  -  (-y,)/y2)1/3c21/c„ 

(B-24) 

s  =  [(C11/C30)2  -  1  ]/ (Cj i/C20)2  - 1 

(B-25) 

s  =  (C20/C2i)2[(C2i/Ch)2  -  1]/[(C20/Cn)2  -  1] 

(B-26) 

and 

t1  =  P 1/ P  2  - 

(B-27) 
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Since  there  is  only  one  value  each  for  these  canonical  parameters,  we  have  deleted 
subscripts  on  the  parameters  to  simplify  the  notation.  Note  that  the  eight-profile  pa¬ 
rameters  reduce  to  four  canonical  parameters. 

The  eigenvalue  equation  is  given  by 


Ai(-  x) 

M  Ai(-  y) 
pAi(-  y) 
0 

where 


Bi(-  x)  0 

/tBi(-y)  -  Ai(-  £21) 

pBi'(-y)  Ai'(-£21) 

0  Ai'(-£2o) 


0 

-  Bi(-  £21) 
Bi'(-?2i) 
Bi'(-  ko) 


(B-28) 


£21  =  P  2(1  -S)x+  p2  S  y 


(B-29) 


and 


£20  =  P  20  ~  S)x+ p2  S  y  (B-30) 

Equation  B-28  follows  from  the  boundary  and  interface  condition  of  eqs.  18  to  24  of 
reference  2  with  appropriate  changes  as  previously  discussed. 

Now,  the  canonical  expressions  are  already  available  here  for  a  two-layer  bounded 
profile  with  no  sound-speed  discontinuity  and  no  treatment  of  density.  Equation  8  may 
be  obtained  from  eq.  B-28  by  the  following  simplifications 


p  =  p 

(B-31) 

n  =  1 

(B-32) 

and 

r* 

Q. 

II 

tN 

(B-33) 

Because  of  this  similarity  we 
to  replace  p2  y  by  £2,,  by  p,  and 

were  able  to  use  the  results  of  section  2B  by  updating 
p  inserting  H  where  appropriate. 

Following  the  treatment  of  section  2B  for  continuous  profile,  we  introduce  various 
sets  of  equations  which  facilitate  the  expansion  of  the  matrix  of  eq.  B-28  and  its  solu¬ 
tion  by  Newton’s  method.  The  counterparts  of  eqs.  12  to  15  become 
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Gi=M  Ai(-  y)Ai'(-  b\)  +  p  Ai'(-y)Ai(-£2i)  (B-34) 

G2  =  /*Bi(-y)Ai'(-£21)  +  P  Bi'(-  y)Ai(-  bi)  (B-35) 

G3  =  /x  Ai(-y)Bi'(-£21)  +  p  Ai'(-  y)Bi(-  bi)  (B-36) 

and 

G4=//Bi(-y)Bi'(-£2I)  +  pBi'(-y)Bi(-C21  .  (B-37) 

The  counterparts  of  eqs.  20  to  23  are  given  by 

G„  =  Bi'(-  bo)  G2  -  Ai'(-  bo)  G4  (B-38) 

G„  =  Bi'(-  bo)  G,  -  Ai'(-  bo)  G3  (B-39) 

G13  =  Bi(-  bo)  G2  -  Ai(—  bo)  G4  (B-40) 

and 

014  =  Bi«2„)  G,  -  Ai(-  ?2„)  Gj  ■  (B-41) 

The  counterparts  of  eqs.  24  to  26  are  given  by 

Gls-Ai(-x)G„  -  Bi(-  x)  Gi2  (B-42) 

G16  =  Ai'(-x)  G„  -  Bi'(-x)  G12  (B-43) 

Gn  =  Ai(-  x)  G13  -  Bi  (-  x)  G14  (B-44) 

The  eigenvalue  equation  for  the  free-rigid  boundary  of  figure  B-l  is  given  by  eq.  28, 
i.e., 

G15  =  0  .  (B-45) 


To  solve  eq.  B-45  by  Newton’s  method  we  need  the  partial  derivatives  with  respect 
to  y  and  with  respect  to  x. 

The  first  step  in  evaluating  these  partial  derivatives  is  to  evaluate  the  partial  deriva¬ 
tives  of  eqs.  B-34  to  B-37.  These  derivatives  are  the  counterparts  of  eqs.  16  to  19  and 
are  given  by 

Gj'  =  -  Ou  +  p  3  S)Ai'(-  y)  Ai'(-  b\)  +  (m  p2  S  bi  +  P  y)  Ai(-  y)  Ai(-  C21)  (B-46) 

G2'  =  -  0u+  p3  S)Bi'(-y)  Ai'(-  bi)  +  (a  p2  S  b\  +  py)  Bi(-y)  Ai(-  bi)  (B-47) 

G3  =  -  (/z+  p3  S)Ai'(-y)  Bi'(-£21)  +  {ji  p2  S  b\  +  py)  Ai(-y)  Bi(- £21)  (B-48) 
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and 


04'  =  -  (M  +  P 3  S)Bi'(-  y)  Bi'(-  $21)  +  (p.  p2  S  £21  +  p  y)  Bi(-  y)  Bi(-  £21)  (B-49) 

The  partial  derivative  with  respect  to  y  is  the  counterpart  of  eq.  30  and  may  be 
written  as 


dGis/dy  =  Gj5'  +  £20p 2  S  Gn  (B-50) 

where  G15'  is  G15  with  G1t  G2,  G3l  and  G4  replaced  respectively  by  G/,  G2',  G3',  and 
G/. 

The  partial  derivative  with  respect  to  x  involves  the  partial  derivatives  of  G,  and  G4 


with  respect  to  x.  Hence,  we  determine  that 

dGj/dx  =  -  p2(1  - S)K2i  Ai (—  y)  Ai(-  £21)  +  Ai'(- y)  Ai'(-  £2 1)]  (B-51) 

<5G2/<5x  =  -  p2(1  -  S)[/^21  Bi(-  y)  Ai(-  £21)  +  Bi'(-  y)  Ai'(-  C21)]  (B*52) 

<5G3/(3x  =  -  p2(1  - S)K21  Ai(- y)  Bi(-  C21)  +  Ai'(- y)  Bi'(-  £2,)]  (B‘53) 

(5G4/ <5x  =  -  P(1  -  S )\m^  Bi(-y)  Bi(-  ?21)  +  Bi'(-  y)  Bi'(-  ?21)]  (B-54) 

The  partial  derivative  with  respect  to  x  may  be  written  as 

(5G15/<5y  =  -  G16  +  ?2o  P 2  (1  -  S)  G17  +  G15  (B-55) 


where  G]5  is  G,5  with  G„  G2,  G3,  and  G4,  replaced  respectively  by  <5Gi/<5x, 
dG2/<5x,  <3G3/5x,  and  <5G4/<5x.  We  note  that  G15  does  not  appear  in  eq.  29,  the  coun¬ 
terpart  of  eq.  B-55.  The  reason  is  that  eqs.  12-15  are  independent  of  x  and  thus, 
eqs.  B-51  to  B-54  are  all  zero. 

One  may  solve  eq.  B-45  by  either  selecting  a  fixed  value  of  y  and  iterating  to  deter¬ 
mine  x  by  eq.  B-55  or  by  selecting  a  fixed  value  of  x  and  iterating  to  determine  y  by 
eq.  B-50.  Iteration  on  y  or  x  then  works  best  for  the  cases  A  or  B  modes  respectively, 
as  discussed  for  continuous  sound  speeds. 

Iteration  using  finite  differences  rather  than  derivatives  is  possible.  However,  we  do 
not  recommend  this  approach  as  the  eigenvalue  curves  often  have  rapid  changes  of 
slope.  If  a  comprehensive  analysis  of  many  cases  is  to  be  made,  then  programming 
<X3i5/<5x  and  <3Gi5/<5y  will  pay  off  in  the  long  run. 

We  note  that  Fell  Hosmer  of  NOSC  Code  711  has  developed  a  comprehensive  com¬ 
puter  program  to  determine  eigenvalues  for  the  two-layer  continuous  sound-speed 
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profile.  This  program  could  readily  be  modified  to  accommodate  the  discontinuous 
sound-speed  profile  including  density  effects.  For  example,  eqs.  B-38  to  B-44  are  iden¬ 
tical  to  those  of  the  continuous  profile.  The  only  significant  modification  is  to  program 
eqs.  B-51  to  B-55  which  do  not  appear  in  the  continuous  program. 

We  note  also  that  the  forms  to  be  used  for  other  boundary  conditions  appear  in  the 
forms  for  the  free-rigid  case.  For  example,  G17  =  0  and  G16  =  0  give  the  eigenvalue 
equation  for  the  free-free  and  rigid-rigid  boundaries  and  G15  appears  in  some  of  the 
partial  derivatives  of  G17  =  0  and  G,e  =  0.  We  will  not  give  the  free-free  and  rigid-rigid 
forms  for  the  discontinuous  profile,  because  these  conditions  are  of  an  academic  inter¬ 
est.  We  note;  however,  that  the  Hosmer  program  treats  these  boundary  cases  and  could 
easily  be  modified  to  do  so  for  the  discontinuous  profile  model. 

Assume  now  that  the  eigenvalues  have  been  determined.  The  normalized  mode 
depth  functions  are  given  by 

Un  (Z0)  Un  (Z)  =  p  (Z0)  p  (Z)  Fn  (Z0)  F„(  Z)/Dn  .  (B-56) 

Here,  n  is  the  mode  number,  ZQ  is  the  source  depth  and  Z  is  the  receiver  depth. 
Equation  B-56  differs  from  eq.  20  of  reference  B2  in  that  eq.  20  does  not  include  the 
density.  Since  the  density  is  the  same  for  all  modes,  reference  B2  introduces  the  den¬ 
sity  as  a  multiplicative  factor  after  the  modes  have  been  summed.  The  reason  we 
chose  to  include  the  density  in  eq.  B-56  is  so  that  plots  of  the  depth  function  will  be 
continuous  functions  of  depth.  If  the  density  is  excluded  in  eq.  B-56,  then  the  depth 
functions  will  have  a  discontinuous  jump  at  each  layer  interface. 

The  unnormalized  eigenfunctions  are  given  by 

Fi(Z)  =  Bi(- x)  Ai(- x  -  a]  Z)  -  Ai(- x)  Bi(- x  -  aj  Z)  (B-57) 


and 


F2(Z)  =  C2  Ai  [-£21  -  a2(Z-  Z2)j  +  D2  Bi  [-£21  -  a2(Z-Z2)]  (B-58) 

where 

a;  =  2tt2/3  f2//3  (-  £i)1/3/Cj .  (B-59) 

For  source  (receiver)  in  layer  1,  eq.  B-57  and  P,  are  used  in  eq.  B-56.  For  source  (re¬ 
ceiver)  in  layer  2,  eq.  B-58  and  P2  are  used  in  eq.  B-56. 

The  form  of  eq.  B-57  satisfies  the  boundary  condition  that  F,  (0)  =  0.  Our  remain¬ 
ing  task  is  to  give  explicit  computing  forms  for  C2,  D2,  and  Dn.  Equations  77  to  79 
apply  without  change. 
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The  counterparts  of  eqs.  80  and  81  lead  to 


D2  =  -tt  [m  Fi(-y)  Bi'(-£21)  -  PF,«(-y)  Bi(-$2i)] 


(B-60) 


and 


Ej.-irblF^-y)  *'(-&,)  -  PF,.(-y)  Bi(-  fc,)]  (B-61) 

Equations  87  and  89  apply  as  is  and  the  counterparts  of  eqs.  92  and  88  are  given  by 

Db  =  -  PM  1  £2oF22  (-  £20)  (B-62) 


and 


D»1  -  (1  +  f  '■  P  3)F>'2(-  y)  +  (y  +  F  p  &,)  F,2  (-  y) .  M 

Equation  B-63  follows  from  eq.  B-17.  Equations  84  and  85  apply  with  eq.  86  replaced 
by 


J  =  C1|C12C2-2  -1|  \2YX\~'  .  (B-64) 

5.  EXAMPLES  OF  DOUBLE  DUCT  PROFILES 

Figure  B-2  presents  schematics  of  four  bottom-interaction  profiles.  Profiles  1  and  2 
represent  shallow  water  profiles  with  a  negative  gradient  in  the  water  column.  Profiles 
3  and  4  represent  deep  water  profiles  with  a  positive  gradient  in  the  water  column. 
Profiles  1  and  3  represent  a  “fast”  sand  bottom  layer.  Profiles  2  and  4  represent  a 
“slow”  clay  bottom. 

Our  interest  here  is  confined  to  profiles  2  and  4  which  are  double  ducts.  A  full 
treatment  of  these  profiles  represent  a  much  more  complicated  problem  than  the 
examples  of  section  3.  The  canonical  parameters  p  and  S  are  replaced  by  P  and  S, 
which  are  complex  rather  than  real.  Moreover  there  are  two  additional  canonical 
parameters— p  which  is  complex  and  H  which  is  real.  Thus,  the  number  of  canonical 
constants  involved  in  full  treatment  is  seven  as  compared  to  two  for  the  profiles  of  sec¬ 
tion  3. 

We  recommend  a  step  wise  approach  using  simpler  profiles  rather  than  a  full-scale 
treatment  of  all  seven  constants.  We  would  first  isolate  the  effect  of  discontinuous 
sound  speeds  by  letting  the  bottom  parameters  be  real  rather  than  complex  and  setting 
M  =1.  This  reduces  the  problem  to  three  real  constants  P,  S,  and  S.  A  further 
simplification  is  to  delete  the  rigid  surface  at  the  bottom  of  layer  2  and  allow  it  to  be 


Figure  B-2.  Schematics  of  the  four  bottom-interaction  profiles. 

an  infinite  half  space.  This  removes  the  canonical  parameter  S  from  the  problem. 
Moreover,  this  reduces  eq.  B-28  from  a  4x4  to  a  3x3  matrix.  We  originally  believed 
that  case  A  of  figure  1  represented  the  simplest  double-duct  configuration.  However, 
the  afore-mentioned  simplifications  of  profiles  2  or  4  result  in  a  simpler  case  than  fig¬ 
ure  1.  Both  cases  contain  two  canonical  parameters.  However,  the  analysis  of  pro¬ 
files  2  and  4  is  easier  because  it  only  involves  a  3x3  matrix.  Thus,  the  characteristics 
of  these  double  ducts  may  be  simpler  than  presented  here  in  section  HI. 

6.  ITERATIVE  APPROACHES  FOR  COMPLEX  EIGENVALUES 

This  section  discusses  an  approach  to  the  solution  of  the  eigenvalue  equation  when 
the  eigenvalues  are  complex. 
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In  the  discussion  of  eq.  60,  reference  1  noted  that  there  must  be  something  about 
the  mathematical  form  of  the  eigenvalue  equation  that  leads  to  Im  x  =  Im  y.  We  deter¬ 
mined  this  property  from  the  physical  condition  that  the  frequency  be  real.  Upon  fur¬ 
ther  thought  we  believe  that  this  is  not  true  about  the  eigenvalue  equation,  i.e.,  there 
are  solutions  for  which  Im  x  /  Im  y.  However,  the  approach  to  treating  this  problem 
is  to  impose  the  condition  Im  x  =  Im  y  as  a  constraint  on  the  problem. 

We  arrived  at  this  belief  by  considering  the  iteration  process  used  to  solve  the 
eigenvalue  equation.  Suppose,  for  example,  we  have  a  problem  that  has  complex 
eigenvalues.  Suppose  further  we  treat  it  by  the  same  iterative  process  that  we  use  for 
real  eigenvlaues,  i.e.,  we  select  a  value  of  y  (real)  and  use  G  and  6G/<5x  to  iterate  to 
the  complex  value  of  x.  However,  this  solution  will  not  be  valid  from  a  physical  stand¬ 
point  because  the  frequency  would  be  complex.  Similarly,  if  we  arbitrarily  picked  y  as 
a  complex  number,  we  would  find  complex  values  of  x  which  satisfy  the  eigenvalue 
equation  but  they  again  would  lead  to  complex  frequency. 

The  trick  to  resolving  this  problem  is  to  insert  the  constraint  into  the  iteration  proc¬ 
ess.  Whenever  we  get  a  new  estimate  of  Im,  x,  we  always  set  Im  y  to  this  same  esti¬ 
mate  prior  to  evaluating  G  and  dG/dx .  Suppose  that  we  already  have  a  valid  solution 
for  some  complex  x  and  complex  y  and  we  wish  to  determine' the  solution  for  stepping 
to  a  new  value  of  Re  y.  We  estimate  the  new  Re  x,  and  new  Im  x  =  Im  y  based  on  the 
previous  solution. 

We  then  evaluate  G  and  divide  by  <5G/<5x  to  obtain  the  incremented  value  of  Re  x 
and  Im  x  to  add  to  the  old  estimate  and  obtain  a  new  estimate  of  Re  x  and  Im  x.  At 
this  stage,  there  are  two  possible  approaches.  The  simplest  is  to  use  this  new  estimate 
of  Im  x  and  set  Im  y  to  this  estimate.  A  second  and  probably  superior  approach,  is  to 
set  Im  x  =  Im  y  equal  to  the  average  value  of  the  last  two  estimates  of  Im  x.  This 
latter  approach  takes  into  account  that  the  new  estimate  of  Im  x  is  predicated  on  the 
previous  estimate  of  Im  y.  This  newer  estimate  cannot  be  correct  if  we  change  Im  y. 
Our  guess  is  that  the  solution  for  Im  x  (and  Im  y)  lies  between  previous  Im  x  and  the 
estimate  of  Im  x  for  the  current  iteration. 

We  can  anticipate  this  process  will  take  somewhat  longer  to  converge  than  for  real 
eigenvalues  because  we  are  moving  Im  y  at  each  iteration.  However,  we  think  that  Re 
x  and  Re  y  are  the  dominant  terms  of  the  solution  and  that  Im  x  and  Im  y  are  rela¬ 
tively  small  corrections.  These  assumptions  may  be  wrong  in  some  cases  with  the  proc¬ 
ess  failing  to  converge. 

We  are  looking  forward  to  a  computer  implementation  for  solving  canonical  eigen¬ 
values  when  the  solutions  are  complex.  There  are  two  outstanding  questions  to  be 
answered.  The  first  question  is  whether  the  iterative  technique,  which  constrains 
Im  x  =  Im  y,  will  converge  under  usual  circumstances.  The  second  question  is  the 
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nature  of  solutions  for  Im  y  *  I m  x.  For  example,  if  we  force  Im  y  =  0,  will  we  find  a 
full  suite  of  curves  of  real  y  versus  complex  x  that  correspond  to  modal  solutions?  Are 
there  physical  problems  that  may  be  represented  by  these  solutions,  or  are  the  solu¬ 
tions  merely  an  exercise  in  pure  mathematics? 

In  closing,  we  note  that  we  have  illustrated  the  iteration  process  for  a  fixed  y  and 
the  use  of  dG/dx  to  determine  x.  A  similar  process  is  recommended  for  a  fixed  x  and 
the  use  of  <5G/<5y  to  determine  y.  Here  again,  the  trick  is  to  constrain  Im  x  =  Im  y 
during  the  iteration  process  regardless  of  whether  we  are  iterating  on  x  or  on  y. 
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